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A  bstract 

This  thesis  develops  and  evaluates  a  number  of  new  concepts  and  tools  for  the  anal¬ 
ysis  of  signals  using  variable  overlapped  windows  and  orthonormal  bases.  Windowing, 
often  employed  as  a  spectral  estimation  technique,  can  result  in  irreparable  distortions  in 
the  transformed  signal.  By  placing  conditions  on  the  window  and  incorporating  it  into 
the  orthonormal  representation,  any  signal  distortion  resulting  from  the  transformation 
can  be  eliminated  or  cancelled  in  reconstruction.  This  concept  is  critical  to  the  theory 
underpinning  this  thesis. 

As  part  of  this  evaluation,  a  tensor  product  based  general  7V-point  fast  Fourier  trans¬ 
form  algorithm  was  implemented  in  the  DOD  standard  language,  Ada.  The  most  prevalent 
criticism  of  Ada  is  slow  execution  time.  This  code  is  shown  to  be  comparable  in  execution 
time  performance  to  the  corresponding  FORTRAN  code.  Also,  as  part  of  this  thesis,  a 
new  paradigm  is  presented  for  solving  the  finite  length  data  problem  associated  with  filter 
banks  and  lapped  transforms.  This  result  could  have  significant  importance  in  many  Air 
Force  applications,  such  as  processing  images  in  which  the  objects  of  interest  are  near  the 
borders.  In  addition,  a  limited  number  of  experiments  were  performed  with  the  coding  of 
speech.  The  results  indicate  the  lapped  transform  evaluated  in  this  thesis  has  potential  as 
a  low  bit  rate  speech  coder. 
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DIGITAL  SIGNAL  PROCESSING  USING  LAPPED 


TRANSFORMS  WITH  VARIABLE  PARAMETER 
WINDOWS  AND  ORTHONORMAL  BASES 

/.  Introduction 

Digital  signal  processing  is  a  vast  engineering  discipline  primarily  concerned  with 
the  representation  of  signals  as  a  sequence  of  numbers  or  symbols,  and  the  subsequent 
processing  of  these  numbers  toward  some  goal  (25).  One  purpose  for  this  processing  is  to 
transform  a  signal  into  a  form  which  is  in  some  sense  more  effectively  transmitted  as  part 
of  a  communications  system.  Decomposition  is  one  method  of  representing  a  signal  in  this 
more  desirable  form,  while  the  recovery  of  the  signal  from  this  form  is  called  reconstruction. 
Figure  1.1  depicts  such  a  system.  The  actual  form  of  the  decomposition  and  reconstruction 
algorithms  used  in  this  thesis  will  be  outlined  in  the  next  chapter. 


Figure  1.1.  Communications  System 


1.1  Windows  in  Digital  Signal  Processing 

The  Fourier  transform  is  indeed  the  cornerstone  of  modern  signal  processing  and 
involver  representation  of  a  signal  using  an  orthogonal  sinusoidal  basis  set.  In  practice, 
ev  ry  signal  to  be  processed  must  be  of  finite  length,  and  this  length  defines  the  minimum 
frequency  spacing  required  for  adjacent  basis  elements  of  a  Fourier  representation.  When 
the  data  is  transformed,  any  frequencies  not  commensurate  with  one  of  the  basis  frequen¬ 
cies  will  exhibit  non-zero  projections  on  all  basis  frequencies  (14:52).  This  phenomenon 
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is  known  as  spectral  leakage.  If  the  primary  goal  of  a  system  is  harmonic  analysis,  spec¬ 
tral  leakage  is  especially  incommodious  and  can  be  lessened  through  the  use  of  windows 
specifically  designed  for  this  purpose.  Harris  provides  a  thorough  discussion  on  this  sub¬ 
ject  in  (14).  When  windows  are  used  as  a  step  in  decomposing  a  signal  for  eventual 
recovery,  however,  the  properties  and  parameters  of  the  windows  to  be  used  need  to  be 
quite  different  from  those  designed  specifically  to  limit  spectral  leakage.  While  spectral 
leakage  may  be  apparent  in  the  decomposed  signal,  the  primary  concern  of  the  decompo¬ 
sition/reconstruction  system  is  proper  reconstruction  of  the  original  signal.  Toward  this 
end,  Princen  and  Bradley  (27)  introduced  a  technique  whereby  the  window  is  an  integral 
part  of  the  orthonormal  basis.  This  less  traditional  use  of  '..  indows  will  be  shown  to  be  an 
important  concept  toward  the  development  of  this  thesis. 

1.2  Perfect  Reconstruction 

Any  algorithm  which  is  developed  as  part  of  a  decomposition/reconstruction  scheme 
attempts  to  reconstruct  a.  signal  which  is  as  close  as  possible  to  the  original  input.  If 
all  decomposition  coefficients  are  used  in  the  reconstruction,  and  the  reconstructed  signal 
is  an  exact  copy  (with  the  possibility  of  a  delay)  of  the  original  signal,  then  the  system 
satisfies  the  perfect  reconstruction  property  (3-5).  Satisfying  this  property  is  an  important 
requirement  for  any  decomposition/reconstruction  .scheme  (27:11.53),  so  that  coding  dis¬ 
tortion  resulting  from  compression  techniques  have  a  minimal  impact  on  reconstruction  of 
the  original  signal. 

1.3  Research  Goals 

Several  goals  were  established  which  guided  the  research  of  this  thesi  ,\r  of  the 
goals  in  some  way  involve  Mie  Generalized  Lapped  Orthonormal  Transform  (GLOT),  which 
will  I  ?  outlined  in  detail  in  Chapter  2.  The»e  goals  are  outlined  below: 

1.  Implement  the  GLOT  in  software  and  validate  the  theory  outlined  in  Chapter  2. 

including  perfect  reeonstruction  using  variable  size  windows  and  overlap. 
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2.  Investigate  use  of  the  GLOT  in  processing  finite  length  signals  without  the  boundary 
effects  associated  with  traditional  approaches. 

3.  Investigate  use  of  the  GLOT  in  speech  processing.  Specifically,  determine  the  com¬ 
pression  capability  of  the  GLOT  using  various  spectral  compression  techniques. 

The  first  goal  is  necessary  because  this  is  the  first  research  regarding  the  GLOT.  After 
implementing  the  complete  software  program,  perfect  reconstruction  will  be  validated  for 
different  size  windows  and  overlap.  Processing  finite  length  signals  often  results  in  dis¬ 
tortions  at  the  signal  boundaries.  The  second  goal  is  to  propose  a  method  of  achieving 
perfect  reconstruction  for  this  type  of  signal.  Finally,  the  third  goal  is  a  logical  extension 
of  previous  work  invloving  speech  and  lapped  transforms,  most  notably  by  Malvar  (18). 

1.4  Speech  Compression  -  An  Application 

One  area  which  will  be  investigated  is  compression  of  speech  data,  since  something 
other  than  a  simple  analog-to-digital  conversion  is  required  to  make  transmitting  speech 
economical  in  terms  of  bit  rate.  Compressing  speech  requires  a  three  way  trade-off  among 
the  goals  of  preserving  intelligibility  and  quality,  limiting  the  bit  rate,  and  minimizing  com¬ 
putational  complexity  (26:225).  The  impact  of  any  compression  technique  on  these  three 
often  incompatible  goals  must  be  continually  assessed  in  determining  the  most  promising 
compression  algorithm.  Transform  coding  is  one  of  the  major  low  bit  rate  speech  cod¬ 
ing  techniques  being  investigated  by  the  military,  since  channel  bandwidth  is  extremely 
limited  in  many  strategic  and  tactical  systems  (38).  This  thesis  will  include  preliminary 
experiments  u.sing  the  GLOT  that  will  be  useful  in  describing  the  design  space  for  a  low 
bit  rate  speech  encoder. 

1.5  Outline 

This  chapter  provides  an  introduction  to  this  thesis,  describing  the  goals  of  this  thesis 
and  particular  field  of  research  within  the  vast  discipline  of  digital  signal  processing.  The 
background  and  theory  which  form  the  basis  for  this  thesis,  up  to  and  including  the  GLOT, 
is  outlined  in  Chapter  2.  Chapter  3  provides  a  complete  derivation  and  analysis  of  an  Ada 
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based  fast  Fourier  transform  algorithm,  the  most  significant  coding  effort  of  this  thesis. 
Chapter  4  provides  the  link  between  the  theory  of  Chapter  2  and  the  application  of  this 
theory  to  various  signal  processing  and  communications  problems,  including  finite  length 
processing,  and  speech  processing  and  compression.  Finally,  a  summary  of  key  results,  and 
a  recommendation  for  follow-on  research  areas,  are  outlined  in  Chapter  5. 
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II.  Background 


In  (33),  Suter  and  Oxley  provide  the  derivation  of  a  new  method  for  the  decompo¬ 
sition  and  reconstruction  of  signals  based  on  variable  parameter  windows  and  weighted 
orthonormal  basis  functions.  The  aspect  of  this  derivation  which  is  different  from  the  con¬ 
ventional  methods  of  transform  coding  is  that  the  orthonormal  basis  used  to  represent  the 
signal  of  interest  is  actually  the  product  of  a  window  and  some  type  of  orthonormal  function 
(in  this  application,  a  sinusoidal  function).  As  described  in  Section  1.1,  this  technique  was 
developed  to  overcome  some  of  the  limitations  of  finite  length  processing.  Malvar  (18)  (21) 
has  investigated  such  a  method,  called  the  Lapped  Orthogonal  Transform  (LOT),  and 
adapted  this  approach  to  speech  signals  to  overcome  what  he  has  termed  “blocking  ef¬ 
fects.”  Combining  the  window  and  function  to  form  an  orthonormal  basis  leads  to  an 
overlap  of  adjacent  intervals,  and  Malvar  limited  his  analysis  to  constant  length  windows 
with  a  50//:%  overlap,  for  positive  integer  k  (19).  He  was  able  to  eliminate  the  “extrane¬ 
ous  tones”  in  speech  resulting  from  the  conventional  windowing  scheme  (21:553).  Coifman 
and  Meyer  (8)  generalized  Malvar’s  analysis  to  include  variable  size  windows.  Suter  and 
Oxley’s  approach  is  based  in  part  on  these  developments,  and  further  generalizes  previous 
work  to  include  variable  bases  (including  non-sinusoidal  bases),  variable  length  windows, 
variable  overlap  between  adjacent  windows,  and  variable  window  amplitude.  They  also 
provide  a  computationally  efficient  algorithm  (33).  The  transformation  outlined  by  Suter 
and  Oxley  will  be  mathematically  defined  in  the  next  section. 

Although  not  an  integral  part  of  this  thesis,  it  is  important  to  note  that  there  are 
intimate  relationships  between  multirate  filter  banks,  lapped  orthogonal  transforms,  and 
wavelets  (2)  (9)  (20)  (35). 

The  initial  implementation,  which  forms  the  basis  for  this  thesis,  is  restricted  to  a 
sinusoidal  basis  with  a  constant  amplitude  window  of  variable  size  and  a  variable  per¬ 
cent  overlap  with  adjacent  windows.  These  initial  assumptions  will  be  used  to  tailor  the 
derivation  of  the  transformation  outlined  in  (33)  to  the  work  encompassed  herein. 
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2.1  The  Generalized  Lapped  Orthonormal  Tran'>form 

2.1.1  Mathematical  Definition  The  following  is  a  formal  definition  of  the  lapped 
transform  used  in  this  thesis. 

Let  N  represent  the  natural  numbers  and  Z  represent  the  integers. 

Let  {ojU  G  Z}  be  a  sequence  such  that  aj  oo  as  j  -*  oo  and  aj  — oo  as  j  — >  — oo. 

Let  {ej\j  G  Z}  be  a  sequence  such  that  Cj  >  0  and  cj  +  <  Oj+i  —  aj. 

Let  {fj,k\k  G  N}  be  an  orthonormal  set  on  [aj,aj+i]  with  respect  to  weight  Pjix). 

Let  Wj{x)  be  a  window. 

Then,  define  the  Generalized  Lapped  Orthonormal  Transform  (GLOT)  T  :  ^  P{Z  X  N). 

[Ts]j,i  =  /  s{x)Uj^k{x)dx,  (2.1) 

J  a  j-cj 

where  Uj^k{x)  =  Wj{x)Jj^k{x)^Pj{x). 

The  form  of  Wj{x),  fj^ic{x),  and  pfix)  used  in  this  thesis  will  be  given  when  describing 
the  signal  decomposition  in  the  following  sections. 

2.1.2  Signal  Decomposition  The  decomposition  portion  of  the  GLOT  involves,  at 
the  highest  level,  three  main  steps  in  the  order  outlined  below: 

1.  Multiplication  of  the  Signal  by  a  Window 

2.  Folding  of  the  Data 

3.  Performing  a  Fourier  Transform 

Each  of  these  steps  will  be  outlined  followed  by  a  summary  of  the  full  forward  decomposi¬ 
tion.  First,  however,  some  notation  must  be  provided.  Each  interval  or  block  of  data  to  be 
processed  will  span  the  closed  interval  [aj,«j+i]  as  shown  in  Figure  2.1.  An  orthonormal 
basis,  {fj  ic{x)\k  G  N},  will  be  defined  on  this  interval.  This  orthonormal  basis  must  satisfy 
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Oj  +(j 


Figure  2.1.  Partitioned  Axis 


the  orthonormality  property  (33),  given  in  Equation  2.2: 


/  fj.ki^)fjA^)Pj{^)dx  =  h.i- 

J  a. 


(2.2) 


where  j  indexes  the  current  interval,  /,,j:(x)  specifies  the  kth  basis  function  within  the  jth 
interval,  and  Pj{x)  specifies  a  weight  function  over  interval  j.  The  weight  function  (and  its 
extension,  Pj(x))  for  a  sinusoidal  basis  function  is  unity  for  all  x.  Because  only  a  sinusoidal 
basis  was  implemented  in  this  thesis,  the  weight  function  may  be  omitted  from  further 
equations. 

The  orthonormal  basis  used  in  this  research  was  the  generalized  sinusoidal  basis  of 
Coifman  and  Meyer  (8)  defined  by  /j,*(x)  on  the  interval  [oj, Oj+i): 


(2.3) 


The  first  three  sinusoidal  basis  functions  for  k  —  0,1,2  over  a  normalized  interval  are 
depicted  in  Figure  2.2.  Next,  the  odd  and  even  extensions  of  /;,i(x)  are  defined  on 
(aj  -  (j^aj)  and  (oj+i, +  fj+i ),  respectively,  resulting  in  a  new  function  fj  ic{x),  defined 
by: 


0  ,  -OO  <  X  <  Qj  -  tj 

sin  +  1/2)  ’  Oj  -  f,  <  a:  <  n,+i  +  fj  +  i  (2.4) 


+  f;  +  l  <  X  <  OC. 
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Figure  2.2.  Sinusoidal  Basis  Functions:  A:  =  0,1,2  over  Normalized  Interval 

To  guarantee  that  the  interval  {aj,aj  +  Cj)  does  not  overlap  with  (oj+i  -  €j+i,aj+i),  it 
is  required  that  Cj  +  Cj+i  <  -  a^.  That  is,  if  cj  =  Cj+i,  then  the  maximum  value  of 

€j  is  50%  of  the  length  of  the  window.  A  general  proof  that  the  product  of  the  window 
and  function  (where  both  satisfy  given  criteria)  form  an  orthonormal  basis  may  be  found 
in  (33). 

The  goal  of  this  algorithm  is  to  derive  a  set  of  coefficients  which  represent  the  data 
(coefficients  which  can  then  be  used  to  reconstruct  the  original  signal).  In  other  words,  the 
signal  of  interest  will  be  expanded  in  terms  of  an  orthonormaJ  basis,  {uj  ic\j  G  Z,fc  e  N). 
The  orthonormal  basis  has  been  derived  as  the  product  of  a  window,  Wj,  and  an  expanded 
orthonormal  basis  function,  fj  i,{x).  The  signal  expansion  in  terms  of  the  resulting  or¬ 
thonormal  basis  (33)  is  given  by; 

H  (2-5) 

jezfceN 
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The  forward  transformation  yields  coefficients  aj  k,  which  are  defined  specifically  by: 

rOj+i+tj+i 

O;,*  =  /  s{^)wj{3:)fj,k{x)dx.  (2.6) 

where  s(x)  is  the  data  to  be  represented,  and  /j,t(a;)  is  as  defined  previously  in  Equation  2.4. 
The  window  function,  Wj{x),  will  be  defined  in  the  following  section. 

2. 1.2.1  Multiplication  of  the  Signal  by  a  Window  To  uniquely  specify  the 
GLOT,  a  window  must  be  defined.  The  window  must  have  certain  properties  to  ensure 
that  the  product  of  this  window  and  the  expanded  orthonormal  basis  function  also  form 
an  orthonormal  basis  set.  Properties  first  defined  by  Coifman  and  Meyer  (8)  include: 


(a) 

Wj(x)  =  1 

for 

X  G  (  U j  T  6j  ,  Qj 1 

W 

0 

II 

for 

X  ^  (flj  —  (j,aj^i  +  G+i) 

(c) 

Wj(aj-(T)  =  wj.iiaj+a) 

for 

a  e 

(d) 

wj(x)  +  wj_i(x)  =  1 

for 

X  G  [oj  —  f  j ,  flj  +  ej  ] 

One  window  satisfying  these  properties  (see  Suter  and  Oxley  (33))  is  given  by: 


0 


-00  <  X  <  Oj  —  (j 


sin  (577  {x  -  [aj  -  Cj]}) 


wj(x)  = 


,  ttj  -  Cj  <  X  <  Oj  +  Cj 

1  ,  Oj  -h  (j  <  X  <  Oj^i  — 

(47^  -  K+i  -  9+1]})  ’  Uj+i  -  fj+i  <  X  <  aj,+i +ej+i 

0  ,  aj+i  +  <  I  <  00. 


(2.8) 


This  sinusoidal  window  was  used  for  all  research  in  this  thesis  and  its  asymmetrical 
overlap  with  the  two  adjacent  windows  is  depicted  in  Figure  2.3.  For  these  windows,  the 
polynomial  within  the  braces  is  an  affine  function  of  x.  Coifman  (7)  has  developed  other 
windows  where  the  polynomial  within  the  braces  is  a  higher  degree  polynomial  function 
of  X.  Another  window  satisfying  all  properties  in  Equation  2.7  which  has  the  additional 
advantage  of  being  differentiable  a.t  x  =  aj  —  (j  and  x  =  Oj+i  +  is  defined  by  Suter 
and  Oxley  in  (33). 
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Figure  2.3.  Example  Window 


2. 1.2.2  Folding  of  the  Windowed  Data  An  important  step  in  GLOT  (33)  is 
the  folding  (and  unfolding  discussed  as  part  of  the  reconstruction)  of  the  data  to  form  a 
new  vector,  hj,  as  shown  in  Equation  2.9  below: 


s{x)Wj{x)  -  s(2aj  -  x)wj(2aj  -  x) 
hj(x)  =  {  s(i) 

y  s(x)wj{x)  +  -s(2aj+i  -  x)Wj{2aj+i  -  x) 


Oj  <  X  <  Oj  -{■  Cj 

Oj  +  (j  <  X  <  Oj+x  -  tj+i  (2.9) 
®j  +  i  ~  ^  ^  ^  ”;  +  i- 


The  folding  is  a  natural  consequence  of  the  mathematical  proof  for  the  orthonormal  ba¬ 
sis  (33),  and  has  several  important  implications.  Consider  the  interval  (oy  —  Cj,aj  -f  fj), 
shown  in  Figure  2.4.  To  demonstrate  the  folding,  the  value  hj(aj  -  1)  will  be  a  sum  of  the 
original  data  a.t  Oj  —  I  and  aj  +  I  multiplied  by  the  appropriate  window  value  (in  this  case, 
the  window  value  indicated  by  arrows  in  Figure  2.4).  Data  in  the  interval  (Oj  -  tj,aj  +  (j) 
will  be  included  in  the  processing  of  data  in  intervals  j  and  (j  -  1).  Recall  that  the 
windowed  overlapped  orthonormal  basis  is  designed  to  eliminate  edge  effects.  The  theory 
predicts  it  will  be  possible  to  reconstruct  the  data  points  in  this  region  as  a  weighted  sum 
or  difference  of  two  coefficients  from  adjacent  intervals.  The  reconstruction  algorithm  will 
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Figure  2.4.  Folding  of  Data 


be  outlined  in  more  detail  in  the  Section  2.1.3.  Henceforth,  all  subsequent  processing  of 
the  data  in  the  jth  interval  can  be  confined  to  [a^, Cj+i],  even  though  this  processing  is 
affected  by  a  larger  number  of  data  points,  [a^  -  ej,Oj+i  +  e^+i] 


2. 1.2. 3  Fast  Fourier  Transform  A  fast  Fourier  transform  (FFT)  is  employed 
in  both  the  decomposition  and  reconstruction  stages  of  the  GLOT.  The  form  of  the  FFT 
used  is  given  by: 


M,~\ 


;=o 


(2.10) 


where  vector  Gj  i,  is  the  one-dimensional  transform  of  vector  gj  i,  and  Mj  is  the  length 
of  these  vectors.  Normally,  this  form  is  considered  an  inverse  FFT  based  on  the  positive 
argument  of  the  exponential  kernel.  The  inverse  FFT  used  in  both  the  decomposition 
and  reconstruction  is  of  this  same  form,  and  it  will  subsequently  be  refered  to  simply  as 
an  FFT.  Normally,  the  resulting  complex  vector  from  the  FFT  has  significant  real  and 
imaginary  parts.  In  the  decomposition,  only  the  real  portion  of  the  FFT  output  vector 
contains  any  information.  In  the  reconstruction,  the  complex  FFT  output  is  scaled  and 
only  the  imaginary  portion  of  the  scaled  output  vector  contains  any  information. 
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2. 1.2.4  Decomposition  Algorithm  Summary  With  the  major  steps  in  the  de¬ 
composition  algorithm  outlined,  a  summary  of  the  full  forward  transformation  can  be 
provided.  The  theory  (33)  allows  arbitrary  spacing  between  samples.  This  spacing  has 
been  assumed  to  be  unity  so  that  Xjj  =  Uj  -f  /  for  /  =  0, 1, ...,  Nj  in  the  following  algorithm 
summary: 

(1)  Obtain  data  s(xj  i)  for  aj  —  Cj  <  Xjj  <  aj+i  -f  fj+i, 
that  is,  /  =  —Mj, ..,  -1,0, 1, ...,  Nj  -f 


(2)  Multiply  s{xj  i)  by  window  Wj{xji)  for  /  =  -Mj,...,  -1, 0, 1, ...,  Nj  -t-  Mj+i. 


(3)  Fold  the  sequence  by  defining 

s{Xj  i}Wj{xj  i)  -  s{2aj  -  Xjj)wj(2aj  -  Xjj)  ,  aj  <  Xjj  <  Oj  +  (j 
Hj.i  —  *  T  ^  ^j,i  ^  ~  ^;  +  i 

s(Xj^i)wj(xjj)  +  s(2aj+i  -  Xjj)wj(2aj+i  -  xjj)  ,  o^+i  -  <  Xjj  <  aj+i 

for  /  =  0, 1, ...,  Nj. 

(4)  Define  new  array 


/3. 


'j.i 


=  for  /  =  0, 1, . . . ,  iVj  . 


(.5)  Define  the  even  extension  of  (ijj 

01.1 


0J.I  = 


,  0<  l  <  Nj  -  l 
0j.2Nj-i  ,  Nj  <  I  <  2Nj  -  1. 


(6)  Perform  an  FFT  of  length  2Nj  on  Pjj  and  define  Dj  k  by: 


2N,-l 


Djk  =  He 


2N-  ^ 
(=0 


Trkll2N, 
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(7)  Interpret  the  results  of  the  FFT 


O!j,o  —  0 

0^j,k  —  ^j,k-l  +  2^2A^j-Dj  i  . 


The  coefficients  after  decomposition,  will  be  referred  to  as  spectral  coefficients 
or  coefficients  throughout  this  thesis.  In  the  case  of  the  speech  application,  any  further 
processing  to  achieve  compression  of  the  signal  will  be  performed  at  this  stage  in  the 
spectral  domain. 

2.1.3  Signal  Reconstruction  The  reconstruction  algorithm  is  the  mathematical  in¬ 
verse  of  the  decomposition,  although  the  reconstruction  implementation  does  not  mirror 
a  reverse  ordered  decomposition.  The  reconstruction  (33)  involves  three  steps  as  detailed 
below: 


(1)  Define  odd  extension  of  aj  k  by: 


,  0  <  It  <  -  1 

—  Olj,2N,-k-l  >  Af j  <  <  2Nj  —  1 


(2)  Perform  an  FFT  of  length  2Nj  on  Oj  k 


for  A:  =  0, 1, 2, ... ,  2Nj  —  1  and  define  Hj j  by: 


Hj.,= 


2N,-\ 


2Nj 


te 


2Kkl 

"T7T7 


k=0 


(3)  Reconstruct  the  signal  on  the  interval  [o^, a^+i]  at  the  data  points  Xj  i  =  Uj  +  I  for 
/  =  0,1,2,.  ..,A^j  by  using 


Wj{2aj  - 

Hu 


[  Wj(Xjj)Hjj  -  wj{2aj+i  - 


,  Uj  <  Xj^i  <  aj  +  €j 

1  ~  ^ j  + 1 

’  ^j+i  ^  ^j,i  ^ 
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The  values  Sjj  will  approximate  the  signal  evaluated  at  Xjj,  that  is,  Sjj  ss  s{xj  i). 

There  are  a  few  noteworthy  differences  as  compared  to  the  decomposition  algorithm. 
First,  only  the  imaginary  (rather  than  the  real)  portion  of  the  complex  scaled  FFT  result 
is  of  interest.  Second,  the  unfolding  and  division  by  the  window  are  combined  in  one  step. 
Notice  that  each  reconstructed  point  in  a  region  of  overlap  requires  the  contribution  of  two 
different  windows  and  coefficients.  This  is  consistent  with  the  fact  that  spectral  coefficients 
were  derived  after  folding  windowed  data  points  from  adjacent  windows  (see  Equation  2.9 
above). 

If  all  spectral  coefficients  from  the  forward  transformation  are  used  in  the  reconstruc¬ 
tion,  perfect  reconstruction  should  result.  In  this  case,  Sj j  =  s(Xj^i). 

2.2  Finite  Length  Signal  Considerations 

A  common  signal  processing  problem  involves  proper  reconstruction  of  finite  length 
signals.  When  performing  sub-band  analysis  using  multirate  filter  banks,  for  example,  the 
total  amount  of  data  required  in  the  channel  is  greater  than  that  of  the  finite  length  input 
(and  reconstructed  output)  (16:162).  Assumptions  must  be  made  about  the  data  beyond 
the  signal  extent,  which  may  create  errors  in  the  reconstructed  output.  The  most  common 
methods  being  investigated  are  various  signal  extensions  (37),  including  zero  padding, 
periodic  extension,  symmetric  extension,  and  boundary  value  replication.  In  practice, 
all  signals  must  be  of  finite  extent.  This  problem  is  especially  significant,  however,  in 
applications  such  as  image  coding  (1)  (5)  (23),  time-frequency  analysis  (17),  and  matrix 
computations  (3)  (4). 

The  GLOT  is  also  susceptible  to  this  problem,  since  overlapping  is  employed  and 
data  on  [a_,  -  -f-  fj+i]  is  required  to  process  and  reconstruct  data  on  [0^,0^+!]. 

There  is  currently  no  lapped-transform  based  solution  to  processing  a  finite  length  signal 
which  achieves  perfect  reconstruction  without  storing  additional  data  outside  the  interval 
of  interest,  although  some  approaches  based  on  approximations  have  been  presented  (10). 
This  section  will  outline  a  solution  to  the  problem  which  employs  the  GLOT. 
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Figure  2.5.  Periodically  Extended  Interval  [flj,  Oj+i] 

Before  discussing  a  solution  to  the  finite  length  signal  problem,  three  goals  will  be 
outlined.  First,  perfect  reconstruction  of  the  entire  signal  is  required.  This  goal  assumes 
that  the  data  at  the  boundaries  must  be  preserved  as  well  as  the  rest  of  the  data.  Second, 
the  transform  employed  should  be  nonexpansive  (5),  meaning  no  e^tra  data  must  be  stored 
outside  the  interval  [oj, Oj+i],  thereby  conserving  datastorage  costs.  Third,  the  complexity 
of  the  method  should  be  no  more  than  for  the  GLOT  outlined  in  the  previous  sections. 
These  three  goals  will  be  used  to  measure  the  success  of  the  solution  provided  below. 
Although  no  additional  data  points  will  be  required,  it  will  be  convenient  to  initially 
visualize  the  periodic  extension  of  the  window  and  data  before  and  after  the  finite  interval. 
Now,  treat  the  boundary  points  of  the  interval  [cj,  o^+i]  as  midpoints  of  the  transition 
regions  of  the  same  percent  overlap.  As  illustrated  in  Figure  2.5,  there  is  a  correspondence 
between  the  overlap  at  the  two  boundaries.  In  other  words,  the  following  window  segments 
are  identical  and  span  an  interval  where  the  data  is  also  identical:  a  and  e,  b  and  /,  c  and 
g,  d  and  h.  Because  of  this  correspondence,  only  the  window  segments  c,  «?,  e,  and  / 
and  data  from  Oj  to  Oj+i  are  required  to  reconstruct  data  on  [«;  ,  Oj+i].  For  e.xample,  the 
product  of  window  segment  a  and  the  data  on  [oj  —  €j,aj]  can  be  replaced  by  window 
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segment  e  and  data  on  [«j+i  -  in  the  GLOT  equations  given  above.  This 

satisfies  the  second  goal  of  a  nonexpansive  transform.  It  should  be  noted  that  perfect 
reconstru''tion  can  be  achieved  by  storing  data  other  than  the  periodic  extension  outside 
the  interval  of  interest.  One  example  would  be  to  zero-pad  outside  [aj,aj+i].  Perfect 
reconstruction  could  be  achieved,  however,  there  would  be  an  additional  requirement  to 
store  the  zeroes  at  specific  locations.  There  is  nothing  inherent  in  this  method  which  limits 
the  reconstruction  capabilities  of  the  GLOT,  including  the  potential  discontinuity  in  the 
data  at  the  boundary.  Also,  the  complexity  of  this  solution  is  no  greater  than  that  for 
the  GLOT.  All  of  the  same  transformation  steps  are  followed,  but  with  data  on  different 
intervals.  This  solution  should  satisfy  all  three  goals  outlined  above.  Not  only  is  this  idea 
a  fairly  simple  theoretical  concept,  it  will  be  shown  in  Chapter  4  to  be  a  very  practical 
approach  to  solving  the  problem. 

2.3  Independent  Interval  Processing 

The  process  outlined  above  solves  the  problem  of  processing  the  first  and  last  seg¬ 
ments  of  a  finite  length  record  of  data.  Tin  ..oncept  could  be  expanded  and  applied  to 
each  interval  of  data,  rather  than  just  at  the  signal  boundaries.  If  the  window  and  data  in 
the  jth  interval  were  periodically  extended  as  shown  in  F'igure  2.5,  then  processing  the  jth 
interval  would  be  independent  of  the  data  in  either  adjacent  window.  Such  a  system  would 
provide  a  decomposition/reconstruction  algorithm  achieving  perfect  reconstruction  on  an 
arbitrary  length  data,  without  requiring  any  overlap  between  adjacent  intervals.  This  type 
of  nonexpansive  signal  processing  has  not  previously  been  defined  for  any  lapped  trans¬ 
form.  Chapter  4  will  provide  the  results  of  using  the  GLOT  for  tliis  type  of  independent 
interval  processing. 
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III.  Development  and  Analysis  of  an  Ada  Based  Fast  Fourier  Transform 


The  most  obvious  and  important  initial  coding  task  was  development  of  a  fast  Fourier 
transform  (FFT)  routine  for  a  general  number  of  input  points,  in  accordance  with  step  6, 
Section  2. 1.2.4,  and  step  2,  Section  2.1.3.  This  routine  was  seen  as  the  most  formidable  in 
terms  of  software  implementation,  and  was  consequently  the  highest  coding  priority.  This 
chapter  outlines  the  implementation  of  the  FFT  used  in  this  thesis,  and  is  based  on  the 
work  by  Raduenz,  et  al.  in  (29). 

3. 1  Choice  of  a  Software  Language 

Because  the  choice  of  a  language  in  thesis  software  development  is  significant,  a  few 
remarks  on  this  choice  are  in  order.  All  software  developed  during  the  course  of  this  thesis 
was  programmed  in  the  high-level  language  Ada.  Ada  was  designed  and  developed  by 
the  Department  of  Defense  in  the  late  1970’s  and  early  1980’s  to  ensure  sound  software 
engineering  concepts  were  employed  in  the  development  of  military  systems.  The  Ada 
advantages  of  modularity,  readability,  testability,  and  reliability  were  paramount  to  the 
success  of  the  coding  phase  of  this  thesis.  Further,  follow-on  students  or  any  interested 
programmer  may  use  the  packages  developed  for  this  thesis,  even  by  calling  them  from 
another  language  such  as  C  or  FORTRAN. 

3.2  Background 

Using  tensor  analysis,  Ferguson  (12)  presented  an  elegant  derivation  of  Classman’s 
(13)  general  N  fast  Fourier  transform  (FFT),  together  with  a  concise  FORTRAN  program 
to  implement  this  FFT.  Moreover,  most  FFT  routines  that  were  developed  after  Ferguson’s 
work  (12)  were  and  continue  to  be  based  on  tensor  analysis  (see,  for  example  Van  Loan 
(.36)).  Nonetheless,  most  FFT  routines  commercially  available  today  operate  on  a  vector 
of  length  N  =  2’^,  where  m  is  an  integer.  The  CLOT,  however,  requires  a  wider  choice  of 
N  for  general  implementation.  One  of  the  practical  advantages  of  Classman’s  routine  is 
that  it  can  be  used  in  digital  signal  processing  applications  for  analysis  of  data  of  arbitrary 
length,  without  the  coding  complexity  of  Singleton’s  case  driven  routine  (30).  The  principal 
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disadvantage  of  Classman’s  routine  is  that  it  requires  an  A^-Vector  working  space.  When 
programmed  in  Ada,  however,  the  implementation  details  of  this  and  other  aspects  of  the 
code  can  be  hidden  in  the  body  of  an  FFT  package  (making  them  transparent  to  the  user), 
as  shown  in  the  next  section. 

Industry  and  laboratories  have  been  somewhat  lethargic  about  commercial  use  of 
Ada,  although  today  the  market  approaches  $1.5  billion  annually  (11).  The  military  has 
also  been  criticized  for  allowing  limitless  waivers  and  generally  impeding  the  transition 
to  Ada.  Many  feel  Ada’s  main  disadvantage  is  slow  execution  time,  thereby  rendering 
it  not  applicable  to  many  engineering  applications.  After  delineating  the  conversion  of 
Ferguson’s  FORTRAN  program  to  Ada,  an  execution  time  comparison  will  demonstrate 
that  Ada  execution  time  is  competitive  with  Ferguson’s  FORTR.4N  FFT  execution  time. 

3.3  Source  Code  Comparison 

Figure  3.1  is  a  listing  of  Ferguson’s  FFT  subroutines.  These  subroutines  were  stored 
and  compiled  as  separate  files.  Figure  3.2  gives  the  analogous  Ada  code  in  a  FFT  package 
specification  and  body.  The  source  code  in  the  two  languages  is  dissimilar  in  a  few  notable 
areas.  First,  the  FORTRAN  program  allows  implicit  assignment  from  a  single  index  array 
to  a  three  index  array  within  the  G/assman  subroutine.  This  assignment  is  handled  by  the 
compiler.  In  Ada,  the  same  assignment  could  be  coded  as  a  triple  nested  loop,  assigning 
each  element  in  Dato^  Vector  to  the  appropriate  location  in  a  new  three  dimensional  vector. 
A  faster  and  more  elegant  assignment  can  be  used,  however,  if  the  appropriate  element  in 
the  single  index  vector  can  be  accessed  based  on  the  values  of  the  three  indices.  It  can 
be  shown  that  some  variable  Data.  l''ector.3/ndex(i,  j,  k),  for  i  G  [1,  .4],j  G  [l,Z?].A’  G  [TC"] 
can  be  represented  as  Data.Vector{{i  -  1)*  B  *C  (j  -  1)  +  C  +  A,-)  where  A,  B  and  C  are 
the  integers  passed  to  Glas.sman.  The  above  expression  can  be  simplified  to  : 

Data.V ect.or(({i  —  1)  *  B  +  (j  —  1))  *  C  +  k).  (3.1) 

The  assignments  shown  in  Figure  3.2  were  derived  by  substituting  the  appropriate  indices 
intoequation  3.1.  When  writing  to  the  vector  Temp,  which  becomes  the  output,  the  indices 
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of  the  three  index  vector  used  in  Figure  3.1  match  the  order  of  the  nested  loops,  and  this 
allows  for  an  Ada  assignment  via  a  simple  counter  variable. 

Second,  the  user  of  the  Ada  version  need  only  pass  a  complex  data  vector  and  boolean 
inverse  operator  to  FFT.  Work  space  is  established  within  the  appropriate  subprocedure, 
and  the  vector  length  N  can  be  determined  using  Ada  array  attributes,  thereby  eliminating 
two  of  the  passed  parameters.  Because  work  space  in  Ada  is  not  defined  at  the  top  level, 
there  is  no  need  to  call  Glassman  using  a  boolean  operator  and  alternating  if  statements 
(passing  either  U  or  Work  as  the  input  vector),  as  is  done  in  the  FORTRAN  routine.  The 
only  reason  to  pass  two  arrays  would  be  to  maintain  integrity  of  the  input  data  while 
passing  data  out  as  a  separate  vector. 

Finally,  the  Ada  code  contains  package  calls  to  Complex.Package,  TypeJPackage,  and 
Math  for  standard  complex  number  manipulations,  global  type  declarations,  and  mathe¬ 
matical  operations  respectively,  as  well  as  a  1/A  scaling  at  the  end  of  subprocedure  FFT 
for  the  inverse  transform.  All  timing  analysis  was  performed  using  the  forward  transform 
so  no  scaling  would  be  invoked  in  either  routine.  In  terms  of  floating  point  operations,  the 
Ada  and  FORTRAN  routines  are  equivalent. 

3.4  Execution  Time  Comparison 

Execution  times  for  the  two  programs  were  compared  using  CPU  time  from  the  Unix 
time  command.  Although  elapsed  CPU  times  are  measured  to  the  50th  second  with  this 
facility,  the  exact  execution  time  for  the  FFTs  is  not  of  great  importance  and  is  highly 
machine  dependent.  The  desired  result  was  a  relative  measure  of  FORTRAN  and  Ada 
execution  time  for  digital  signal  processing  algorithms  such  as  the  one  used  here. 

Runs  were  made  on  a  Sun  SPARCstation  2  (operating  system  version  SunOS  4.1.2), 
with  very  light  additional  load.  The  software  packages  used  in  this  comparison  were  Sun 
FORTRAN  (Enhanced  FORTRAN  77)  Version  3.1.1  and  Verdix  Ada  Version  6.0.3(d).  The 
results  are  given  in  Table  3.1.  Notice  that  executables  were  developed  which  provided  no 
output  or  wrote  to  a  file  to  isolate  the  effects  of  I/O  in  the  comparison.  As  is  shown,  the 
Ada  execution  times  are  comparable  to  the  FORTRAN  execution  times. 
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subroutine  FFT-Sub  (N,  U,  Work,  Inverse) 
integer  N 

complex  U(N),  Work(N) 
logic^kl  Inverse 
integer  A,  B,  C 
logical  Inu 
A  =  1 
B  =  N 
C  =  1 
Inu  =.  true. 

if  (B  .GT.  1)  goto  30 
if  (Inu)  return 
do  20  i  =  1,  N 

U(i)  =  Work  (i) 
continue 
return 

A  =  C  »  A 
do  40  C  =  2,  B 

if  (mod(B,C)  .EQ.  0)  go  to  50 
continue 
B  =  B  /  C 

if  (Inu)  call  Glassman(A,  B,  C,  U,  Work,  Inverse) 

if  (  .Not.  Inu)  call  Glassman(A,  B,  C,  Work,  U,  Inverse) 
Inu  =  .Not.  Inu 
go  to  10 

end 


subroutine  Glassman  (A,  B,  C,  Uin,  Uout,  Inverse) 
integer  A,  B,  C 

complex  Uin(B,C,A),  Uout(B,A,C) 
logical  Inverse 

This  subroutine  is  called  from  FFT.Sub 


Complex  Delta,  Omega,  Sum 
Twopi  =  6.28318530717958 

Angle  =  Twopi  /  lloat(A«C) 

Delta  =  CMPLX(Cos(Angle),  -Sin(Angle)) 

if  (Inverse)  Delta  =  CON.IG  (Delta) 


Omega  =  CMPLX(t. 0,0.0) 
do  40  IC  =  1,  C 

do  30  lA  =  1,  A 

do  20  IB  =  1,  B 

Sum  =  Uin(IB,C,IA) 
do  10  JCR  =  2,  C 

JC  =  C  +  1  -  JCR 

Sum  =  Uin(IB,JC,IA)  +  Omega  •  Sum 
continue 

Uout(  IB,1A,IC)  =  Sum 
continue 

Omega  =  Delta  *  Omega 
continue 
continue 


return 

end 


Figure  3.1.  Fortran  Code  for  FFT  Subroutines 
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with  ComplexJPkg; 
with  Type-Package; 
with  Math; 


use  Complex-Pkg; 
use  Type-Package; 
u&e  Math; 


package  FFT-Pack  is 

procedure  FFT  (  FFT-Data  :  in  out  Contplex-Veclor  ; 

Inverse-Transform  :  in  boolean  ); 

end  FFT-Pack; 


package  body  FFT-Pack  is 

procedure  Glassman  (  A,  B,  C  in  integer; 

Data-Vector  in  out  Complex-Vector  ; 

Inverse.Transform  in  boolean  )  is 


Temp 

Counter,  JC 
Two -Pi 

Del.  Omega,  Sum 

Angle 

C_PlU3_l 


CompleK_Vcctor(  l..A*B*C); 
integer  :=  0; 

constant  float  :=  6.28316530717958; 

Complex; 

float; 

integer  :=:  C  1; 


begin 

Angle  :=  Two-Pi  /  (float(A«C)); 

Del  ;=  Complex-Of((Co8(Angle)),  (.(Sin(Angle)))); 

if  (Inverse-Transform)  then 

Del  :=  Conjugate(Del); 

end  if; 

Omega  :=  Complex-Of(  1 .0,0.0); 

for  IC  in  1..C  loop 

for  lA  in  l.  A  loop 

for  IB  in  1..B  loop 

Sum  ;=  Data-Vector((((!A  •  1)*C  +  (C-l))  •  B)  +  IB); 
for  JCR  in  2..C  loop 

JC  ;=  C-Plus-l  -  JCR;  —  No  need  to  add  C  1  each  time  through  loop 
Sum  Data-Vector((((IA  •  1)*C  +  (JC  •  1^)•B)  -f  IB)  +  (Omega  •  Sum) 
end  loop;  —  JCR 
Counter  :s  Counter  +  1; 

Temp(Counter)  :«  Sum; 
end  loop;  —  IB 
jmega  :a  Del  •  Omega; 
end  loop;  —  lA 
end  loop;  —  IC 

Data-Vector  :=  Temp; - assign  output  back  to  Data-Vector 

end  Glassman; 


procedure  FFT  (  FFT-Data  :  in  out  Complex-Vector; 

Inverse.Transform  :  in  boolean  )  is 

A  integer  :s  1; 

B  *eger  :a  FFT-Data'length; 

C  ii.  ■!ger  ;a  1; 

begin  —  FFT 

while  (B  >  1)  loop  —  define  the  integers  A,  B,  and  C 

A  ;=  C  •  A;  —  such  that  A*B*C  =  FFT-Data’length 

C  :=  2; 

while  (B  mod  C)  /=  0  loop 
C  :=  C  +  1; 
end  loop; 

B  :=  B/C;  -  -  B  =  1  causes  exit  from  while  loop 

Glassman  (A,B,C,  FFT-Data.  Inverse.Transform); 
end  loop; 

if  Inverse.Transform  then  —  optional  l/N  scaling  for  inverse  transform  only 
for  i  in  FFT-Dala'range  loop 

FFT-Data(i)  ;=  FFT-Data(i)  /  float(FFT J^ata’length); 
end  loop; 
end  if; 
end  FFT; 
end  FFT-Pack; 


Figure  3.2.  Ada  Code  for  FFT  Package 


No  Output 

Write  to  File 

Average  CPU  Time  in  Seconds 

Points 

Ada 

FORTRAN 

Ada 

FORTRAN 

500 

0.1 

0.1 

0.1 

0.3 

0.2 

0.3 

0.3 

0.6 

KnfiM 

0.4 

0.6 

0.8 

1.3 

3000 

0.7 

0.9 

1.2 

1.9 

4000 

0.9 

1.2 

1.7 

2.6 

5000 

1.2 

1.5 

2.2 

3.3 

6000 

1.5 

1.8 

2.7 

3.8 

7000 

1.8 

2.2 

3.2 

4.6 

8000 

2.0 

2.5 

3.6 

5.2 

9000 

2.3 

2.8 

4.1 

5.8 

10000 

2.6 

3.2 

4.6 

6.5 

11000 

3.2 

3.8 

5.5 

7.4 

12000 

3.1 

3.8 

5.6 

7.8 

13000 

4.0 

4.7 

6.7 

9.0 

14000 

3.9 

4.7 

6.9 

9.3 

15000 

4.0 

4.9 

7.2 

9.9 

Table  3.1.  Ada  vs  Fortran  Execution  Time  Comparison 


Although  run  times  were  roughly  the  same,  the  Ada  compilation  time  was  noticeably 
longer  than  that  for  FORTRAN.  One  reason  for  this  extended  compilation  time  is  that 
Ada  was  developed  as  a  strongly  typed  language,  and  performs  numerous  compilation 
checks  to  minimize  run  time  errors.  Ada  also  performs  run  time  checks  such  as  numeric 
range  checking,  which  can  be  disabled  using  the  -5  command.  Because  Sun  FORTRAN 
3.1.1  does  not  perform  this  level  of  run  time  checking,  the  final  Ada  executable  used  in  this 
comparison  was  compiled  with  run  time  checks  disabled.  A  significant  speed  increase  could 
be  realized  in  the  FORTRAN  program  by  compiling  with  the  -fast  or  -fnonstd  commands, 
however,  this  compilation  results  in  floating  point  outputs  which  do  not  conform  to  IEEE 
standards.  Thus,  neither  the  -fast  FORTRAN  nor  the  -O  Ada  optimization  options  were 
used  during  compilation.  It  is  important  to  note  that  the  accuracy  of  both  programs 
outputs  was  comparable,  since  the  outputs  agreed  down  to  the  roundoff  of  the  machine. 
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3.5  Summary 

The  FFT  is  the  most  involved  algorithm  in  the  GLOT.  Coding  this  algorithm  effec¬ 
tively  in  Ada  provided  confidence  that  the  full  transform  could  be  implemented  in  Ada. 
The  execution  time  comparison  indicates  that  Ada  is  competitive  with  FORTRAN  and  is 
a  viable  language  for  other  digital  signal  processing  applications  as  well. 


IV.  Methodology 


Programming  a  general  ^-point  FFT  routine  was  seen  as  the  first  priority  software 
coding  task,  as  discussed  in  Chapter  3.  Before  a  full  software  implementation  could  be 
achieved,  however,  a  prototype  system  was  developed  to  test  the  theory  of  the  GLOT. 
This  chapter  provides  the  link  between  theory  and  application,  ranging  from  initial  imple¬ 
mentation  of  a  prototype  system  to  application  of  the  GLOT  to  current  signal  processing 
problems.  Personal  Computer  -  Digital  Signal  Processing  (PC-DSP),  discussed  below,  wcis 
the  first  method  used  to  test  the  GLOT.  The  program  was  used  as  a  fast  prototype,  and 
aided  in  development  and  testing  of  the  mathematical  theory.  The  next  test  phase  involved 
a  computer  program  implementing  the  GLOT.  The  development  of  this  program  also  relied 
heavily  on  the  previous  PC-DSP  results.  Some  variations  of  the  GLOT  were  implemented 
and  tested,  each  working  toward  a  different  goal.  Finally,  the  GLOT  was  applied  to  digital 
speech  processing.  The  main  focus  of  this  segment  of  the  research  was  to  determine  the 
minimum  number  of  coefficients  needed  to  reconstruct  a  quality  reproduction  of  the  input 
signal,  while  maintaining  an  acceptable  computational  penalty.  The  later  portion  of  this 
chapter  outlines  the  various  techniques  employed  toward  this  end. 

4.1  A  Prototype  System  -  PC-DSP 

After  carefully  delineating  the  theory  behind  the  GLOT,  the  first  step  in  the  evo¬ 
lutionary  process  of  validating  the  theory  was  to  develop  a  prototype  system.  The  tool 
chosen  for  this  purpose  was  Personal  Computer  -  Digital  Signal  Processing  (PC-DSP), 
an  interactive  menu-driven  software  package  for  use  in  common  digital  signal  processing 
tasks  (24).  This  program  performs  many  of  the  functions  required  to  test  the  GI  DT,  such 
as  fast  Fourier  transforms  (FFTs),  point-by-point  vector  multiplications  and  additions, 
vector  scaling,  and  the  design  and  analysis  of  digital  filters.  Limitations  which  constrained 
the  program  included  the  inability  to  loop  or  perform  an  FFT  on  data  of  length  other  than 
2"’,  where  m  is  an  integer.  Another  disadvantage  of  the  program  was  that  each  step  had  to 
be  manually  performed  for  each  window  of  each  data  .set.  Some  automation  was  available 
by  using  PC-DSP  macro  programs,  where  multiple  PC-DSP  commands  could  be  loaded 
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and  executed  from  a  file.  These  macro  programs  also  had  limitations,  mainly  involving 
insufficient  memory  for  this  application.  In  spite  of  these  limitations,  prototyping  using 
PC-DSP  provided  the  ability  to  determine  the  feasibility  of  the  GLOT,  and  also  limited 
the  time  required  to  code  the  eventual  algorithm  in  a  high-level  language. 

The  signal  shown  in  Figure  4.1  was  chosen  to  test  the  analysis  and  synthesis  al¬ 
gorithms.  It  consists  of  a  sinusoid,  a  decaying  exponential,  a  growing  exponential,  and 
another  slowly  varying  sinusoid  in  the  first  through  fourth  partitions  defined  by  the  fol¬ 
lowing  equations: 


sin(0.041x) 


s{x)  =  i 


0.862ej;p[-0.006(a:-  512)] 
0.0402ea:p[0.006(x-  1024)] 
0.862-1-  1.5sin[0.03(x-  1536)] 


,  0<x<511 
,  512  <  a:  <  1023 
,  1024  <x<  1535 

,  1536  <  X  <  2048 


(4.1) 


The  window  boundaries  were  chosen  to  be  512  and  1536,  resulting  in  three  analysis  intervals 
of  512,  1024,  and  512  points  respectively.  The  window  sizes  were  specifically  chosen  to  be 
a  power  of  two  so  that  no  zero-padding  would  be  required  to  perform  the  subsequent  FFT. 
Such  padding  would  have  created  interpolation  errors.  In  an  effort  to  gain  some  insight 
into  the  effects  of  amount  of  window  overlap  on  the  reconstruction  of  the  signal,  the  overlap 
at  data  sample  512  was  51  points  (5%  of  the  larger  window)  and  the  overlap  at  1536  was 
102  points  (10%  of  the  larger  window). 

The  spectral  coefficients  obtained  after  performing  the  forward  GLOT  are  given  in 
Figures  4.2  through  4.4.  Notice  that  only  the  first  thirty  coefficients  of  each  window  are 
shown.  Although  the  transform  yields  as  many  coefficients  as  there  are  data  points  in  the 
window,  the  coefficients  beyond  thirty  for  this  signal  were  essentially  zero.  Only  the  first 
thirty  coefficients  were  plotted  to  produce  the  coefficients  of  interest. 

Intuitively,  the  spectrum  plots  are  what  would  be  expected  for  each  analysis  interval. 
The  GLOT  can  be  roughly  equated  to  a  Fourier  transform,  as  both  are  based  on  represen¬ 
tation  of  a  signal  using  sinusoidal  basis  functions.  Figure  4.2  has  a  strong  6th  coefficient 
which  results  from  the  sinusoidal  basis  corresponding  to  the  frequency  of  the  sinusoid  in 
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Figure  4.1.  PC-DSP  Input  Signal 

the  first  window.  Because  the  frequency  of  this  sinusoid  and  the  nearest  basis  function 
frequency  do  not  match  e,xactly,  other  low  frequency  bases  are  also  noticeable  but  of  a 
lesser  magnitude  (see  Section  1.1).  Another  source  of  noticeable  low  frequency  spectral 
components  in  the  first  window  is  the  overlapping  window.  A  small  portion  of  the  low 
frequency  exponential  is  characterized  by  the  coefficients  in  the  first  window. 

The  signal  in  the  middle  segment  of  Figure  4.1  is  slowly  varying,  and  the  correspond¬ 
ing  spectrum  ( Figure  4.3)  has  many  more  significant  coefficients  than  Figure  4.2.  In  general 
more  spectral  coefficients  are  required  to  represent  a  slowly  varying  time  signal  which  is 
non-sinusoidal.  By  analogy,  the  Fourier  transform  of  the  unit  gate  or  rect  function  is  a 
sine  function  (31:47,83),  a  spectral  representation  which  also  requires  a  significant  number 
of  damped  oscillatory  components. 

As  in  Figure  4.2,  Figure  4.4  shows  a  dominant  coefficient  (4th)  corresponding  to 
the  frequency  of  the  signal  in  the  third  interval.  Additionally,  the  large  0th  coefficient 
corresponds  to  the  DC  offset  in  the  interval. 

After  applying  the  inverse  transform,  the  output  shown  in  f'igure  4.5  wa.s  achieved. 
The  output  signal  and  input  signal  are  essentially  identical,  with  an  rms  error  of  2.47  x  10“'’. 
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Figure  4.2.  First  30  Spectral  Coefficients  of  the  First  Window 
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Figure  4.3.  First  30  Spectral  Coefficients  of  the  Second  Window 


4-4 


The  use  of  PC-DSP  was  invaluable  to  the  development  of  this  thesis.  It  provided 
a  proof  of  concept  indicating  the  developed  theory  of  Chapter  2  was  correct  and  worthy 
of  further  research  without  fully  coding  the  algorithm  in  a  high-level  language.  Slight 
modifications  and  corrections  to  the  theory  were  also  facilitated,  i.e.  making  inequalities 
in  the  derivation  consistent  for  the  “limiting”  one  point  overlap  case.  Finally,  a  more 
retrospective  benefit  was  that  during  coding  of  the  algorithm,  this  same  example  was 
used.  Having  saved  the  plots  after  each  step  in  the  forward  and  inverse  transform,  the 
expected  result  was  known  and  each  step  could  be  validated  before  proceeding  to  the  next 
stage.  Attempting  to  debug  a  program  of  this  magnitude  without  knowing  intermediate 
results  would  have  been  rather  arduous. 

J,.2  Fast  Fourier  Transform 

After  developing  the  prototype  system  and  validating  the  theory  behind  the  GLOT, 
actual  code  was  developed  which  could  provide  faster,  more  general  results.  The  most 
important  initial  coding  task  was  development  of  an  FFT  routine,  and  a  full  discussion  of 
the  Ada  based  FFT  developed  for  this  thesis  is  detailed  in  Chapter  3. 

4.3  System  Software  Development 

With  a  thoroughly  tested  FFT  routine  available,  the  software  development  unique  to 
the  GLOT  could  be  initiated.  One  of  the  underlying  goals  during  all  software  development 
was  to  ensure  the  code  could  be  easily  understood  and  modified.  The  use  of  Ada  as 
a  programming  language  aided  in  achieving  this  goal.  Also,  all  algorithmic  coding  was 
made  as  general  as  was  practical.  This  thesis  is  the  incipiency  of  research  regarding  the 
GLOT,  and  it  is  essential  that  follow-on  students  understand  the  details  of  the  software 
and  can  use  it  as  a  baseline  for  their  subsequent  work.  Therefore,  this  section  will  briefly 
outline  each  package  or  subprocedure  of  the  basic  software  program  and  how  the  package 
or  subprocedure  relates  to  the  complete  program.  Further,  each  portion  of  the  code  will 
be  referenced  to  the  appropriate  section  or  equation  in  Chapter  2.  A  complete  listing  of 
all  source  code  (with  similar  referencing)  is  provided  in  Appendix  A. 
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4.3.1  Main  Program  Subprocedures 


4. 3. 1.1  Get-File-Data  Procedure  Get-File-Data  prompts  the  user  for  the  name 
of  the  data  file  which  contains  all  the  information  necessary  to  run  the  program.  This  input 
file  must  include  the  following  entries  : 

1.  The  file  containing  the  data  set  to  be  processed 

2.  The  desired  filename  for  the  derivative  data 

.1.  The  desired  filename  for  the  spectral  coefficient  data 

4.  The  desired  filename  for  the  output  data 

.'i.  The  desired  filename  for  the  error  data 

6.  The  desired  filename  for  the  information  pertaining  to  this  run 

7.  The  number  of  points  to  be  processed 

8.  The  number  of  data  partitions 

9.  The  partition  point  and  overlap  size  for  each  partition 

Three  different  stages  in  the  transformation  were  of  interest,  including  the  derivative  of  the 
spectral  coefficients  (used  mainly  in  initial  debugging),  the  spectral  coefficients,  and  the 
output  (or  reconstructed  signal).  This  data  was  stored  in  the  user  defined  files  specified  in 
entries  2-4  above.  The  point-by-point  error  between  the  output  and  input  is  captured  in 
the  file  specified  by  entry  5.  General  information  pertaining  to  the  current  data  run  and 
any  compression  or  filtering  employed  can  be  written  to  the  information  file  specified  by 
entry  6.  After  the  number  of  partitions  is  entered,  the  partition  point  and  number  of  points 
overlap  is  loaded  for  each  partition.  Also  in  this  subprocedure,  the  first  window  of  input 
data  is  copied  to  the  end  of  the  file.  This  allows  so  the  data  to  be  processed  as  a  periodic 
waveform  in  an  effort  to  avoid  boundary  noise  at  the  front  and  back  of  the  reconstructed 
signal.  An  example  input  file  is  given  in  Figure  4.6.  This  particular  input  file  specifies 
that  1000  samples  will  be  extracted  from  file  “input.dat”.  The  results  of  interest  will  be 
written  to  the  appropriate  filename.  The  data  will  be  partitioned  into  4  segments,  with 
overlap  of  10  points  for  the  first  two  intervals,  and  20  points  for  the  last  two  intervals. 
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Figure  4.6.  Example  Input  File 

4-3. 1.2  Multiply. By. Window  Procedure  Multiply. By. Window  creates  a  win¬ 
dow  (defined  by  equation  2.8)  based  on  the  size  of  the  data  to  be  processed  and  the  overlap 
on  either  side  and  multiplies  this  window  by  the  data.  It  also  folds  in  the  appropriate  points 
to  create  a  vector  from  aj  to  aj+i  as  defined  in  equation  2.9. 

4-3. 1.3  Compute.Coefficients  Procedure  Compute.Coefficients  performs  an 
even  extension  of  the  data,  and  an  FFT  as  defined  in  steps  5  and  6  of  Section  2. 1.2.4. 
Note  here  that  the  points  into  the  routine  are  from  Oj  to  Gj+i  or  TV  -|-  1  points,  while  the 
coefficients  out  of  the  routine  are  from  Uj  to  Gj+i  —  1  or  coefficients. 

4.3. 1.4  Reconstruction.FFT  Procedure  ReconstructionJ'FT  is  the  first  step 
in  the  reconstruction  of  the  signal,  involving  an  odd  extension  of  the  data  and  an  FFT  as 
defined  in  steps  1  and  2  of  Section  2.1.3.  The  output  data  is  the  imaginary  output  from 
the  FFT. 


4.3. 1.5  Divide. By.  Windows  Procedure  VFmdouvs  employs  the  win¬ 

dow  created  in  Multiply.By.  Window  to  reconstruct  the  signal  as  defined  in  step  3  of  Sec- 
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tion  2.1.3.  The  output  data  from  Oj+i  —  ej+i  to  cannot  be  reconstructed  during  the 
current  window  processing  loop  because  this  segment  requires  the  processing  of  data  from 
aj  to  Uj  +  Cj  of  the  next  window.  Therefore,  during  any  given  window  processing  loop, 
output  data  from  Oj  —  Cj  to  Oj+i  —  is  calculated  and  written  to  the  output  file. 

4. 3. 1.0  Do-Work  Procedure  Do.Work  performs  the  complete  GLOT  by  call¬ 
ing  the  procedures  discussed  above.  Each  data  segment  is  processed  in  the  main  loop  and 
the  variables  characterizing  the  current  data  segment  are  initialized  as  a  first  step  inside 
this  loop.  Any  processing  of  the  spectral  coefficients  to  obtain  signal  compression  can  be 
implemented  between  procedures  Compute.Coefficients  Reconstruction  J  t-  T 

This  procedure  also  implements  the  solution  to  the  finite  signal  length  problem  out¬ 
lined  in  Section  2.2.  The  first  and  last  data  segments  are  handled  as  special  cases.  The 
reason  for  this  special  processing  is  the  absence  of  data  left  of  the  first  segment  and  right 
of  the  last  segment.  Lack  of  data  in  these  regions  means  non-overlapped  processing  on 
[aj,aj  -t-Cj]  in  the  first  segment,  and  on  [oj+i  —  in  the  last,  leading  to  edge  effects 

in  the  reconstructed  output.  The  problem  is  handled  by  processing  the  first  data  segment 
twice.  Looking  back  to  procedure  Get.Data  discussed  in  Section  4. 3. 1.1,  the  first  window 
of  data  was  written  to  the  end  of  the  input  vector.  When  the  input  data  is  passed  to 
Do-Work  for  processing,  the  first  window  of  data  will  be  processed  once  during  the  first 
processing  cycle  (j  =  0)  and  once  during  the  last  (j  =Partitions’last-l).  The  first  time 
through  the  processing  loop,  there  will  be  perfect  reconstruction  in  (1)  the  region  of  no 
overlap  [aj  -I-  Cj ,  aj^i  -  Cj+i  j  ,  and  (2)  the  overlapped  region  with  the  right  adjacent  interval 
[flj+i  -  fj+i.Oj+i].  The  reconstructed  values  from  [aj,aj  -|-  Cj]  will  be  determined  in  the 
final  loop,  where  the  left  adjacent  window  of  data  is  the  last  input  segment.  These  values 
are  written  to  the  beginning  of  the  output  file  in  the  loop  following  a  call  to  procedure 
Divide. By-Windows  (see  Appendix  A).  The  last  input  data  segment  is  properly  processed 
during  the  second-to-last  processing  loop,  since  the  first  data  segment  has  been  appended 
to  the  right  of  the  end  data.  By  “wrapping”  the  data  in  this  manner,  the  reconstructed  out¬ 
put  for  the  first  and  last  data  segments  is  free  from  edge  effects  caused  by  non-overlapping 
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intervals.  All  other  segments  of  data  are  inherently  immune  to  this  effect  as  discussed  in 
the  theory  of  Chapter  2. 

4.3.2  Package  Type.Package  Package  TypejPackage  contains  all  types  used  in  the 
program,  and  these  types  are  global  to  all  packages  and  subprocedures.  Unconstrained 
arrays  are  used  wherever  possible  to  allow  dynamic  memory  allocation.  Types  Direction, 
Phase-Vector,  and  Phase-Type  were  used  as  part  of  compression  techniques  which  will  be 
discussed  in  subsequent  sections. 

4.3.3  Package  FFT-Pack  Package  FFT-Pack  is  fully  explained  and  documented  in 
Chapter  3. 

4.3.4  Package  Vector-Package  Package  Vector-Package  contains  three  functions 
which  contributed  to  the  readability  of  the  main  program,  and  operated  on  unconstrained 
array  sizes.  Function  Even-Extend  implemented  the  even  extension  defined  in  Section  2. 1.2.4, 
step  5,  and  is  called  from  procedure  Compute-Coefficientsmthin  the  main  program.  Func¬ 
tion  Odd.Extend  implemented  the  odd  extension  defined  in  Section  2.1.3,  step  1,  and  is 
called  from  procedure  Reconstruction J'FT.  Finally,  function  Reverse j\.ssignmentvids  used 
to  assign  window  values  to  the  left  window  segment  from  the  right  window  segment  within 
procedure  Divide-By-Windows.  This  function  exploits  the  fact  that  overlapping  window 
segments  must  be  symmetric  about  aj.  Without  this  function,  an  additional  window  vector 
would  have  had  to  be  passed  out  of  procedure  Multiply -By -Windows  and  into  procedure 
Divide.By  Windows.  This  package  can  be  expanded  to  include  future  vector  manipulation 
requirements,  or  could  be  linked  for  use  by  another  algorithm  requiring  these  functions. 

4.3.5  Package  Complex-Pkg  Package  Complex-Pkg  contmns  basic  complex  number 
manipulations  on  real  and  complex  numbers  and  vectors. 

4.3.6  Independent  Interval  Processing  All  of  the  software  code  referenced  thus  far 
was  based  on  overlapped  adjacent  intervals,  except  for  the  two  intervals  at  the  signal 
boundaries.  These  two  intervals  were  treated  as  special  cases  as  discussed  in  Section  4. 3. 1.6. 
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Applying  the  concept  of  independent  interval  processing  discussed  in  Section  2.3  required 
significant  modifications  to  this  software  program.  All  affected  procedures  are  given  in  Ap¬ 
pendix  C  and  replace  those  given  in  Appendix  A.  The  most  important  distinction  evident 
in  the  changes  is  found  in  Procedure  Do.Work.  Previously,  reconstruction  of  the  signal  in 
the  jth  interval  on  [oj+i  —  was  delayed  until  the  processing  of  the  next  interval. 

This  requirement  is  based  on  the  dependency  between  the  overlapped  region  of  adjacent 
signal  intervals.  Using  the  approach  of  Section  2.3,  however,  each  interval  can  be  fully 
reconstructed  independent  of  any  other  intervals.  The  independent  processing  of  data  can 
be  an  advantage  in  terms  of  real-time  processing,  error  control,  parallel  processing,  and 
decreased  buffer  requirements,  all  of  which  will  be  discussed  in  Chapter  5. 

4-4  Applications 

After  developing  a  software  program  which  implemented  the  GLOT,  the  program 
was  used  to  investigate  various  potential  applications.  The  research  discussed  in  this 
chapter  concerning  speech  involves  various  attempts  to  efficiently  encode  the  decomposed 
spectral  coefficients  in  a  manner  that  allowed  for  a  quality  reconstruction.  Applying  the 
GLOT  to  speech  signals  was  a  natural  extension  of  previous  work  in  this  area  (18:975-976). 
Further,  it  was  assumed  that  speech,  being  sinusoidal  in  nature  (32),  would  lend  itself  to  a 
transform  based  on  sinusoidal  representation.  The  following  sections  detail  initial  testing 
and  evaluation  of  the  algorithm,  including  its  use  in  solving  the  finite  length  data  problem, 
perfectly  reconstructing  independently  processed  intervals,  and  compressing  speech. 

4-5  Initial  Testing 

The  PC-DSP  example  outlined  in  Section  4.1  was  used  throughout  development  of 
the  code  outlined  in  Section  4.3.  After  ensuring  the  signal  defined  in  equation  4.1  could  be 
reconstructed  to  roundoff  error,  the  program  needed  to  be  tested  with  an  actual  signal.  A 
sample  speech  signal  from  a  professional  speech  data  base  (22)  was  used,  which  consisted 
of  the  phrase  “American  newspaper  reviewers  like  to  call  his  place  nihilistic”,  spoken  by 
a  male  speaker.  This  60,080  sample  speech  file  was  sampled  at  16  kHz  and  represented 
using  16  bit  integers.  These  numbers  become  important  when  discussing  data  rates  and 
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Figure  4.7.  Sample  Speech  Signal 

compression  techniques.  Figure  4.7  displays  the  first  10,000  samples  of  the  data  which 
correspond  roughly  to  the  word  “American...”.  The  signal  was  divided  into  equal  sized 
windows  of  512  points  with  a  one  point  overlap.  The  one  point  overlap  was  chosen  to 
determine  an  upper  bound  on  the  amount  of  edge  effects  as  a  function  of  window  overlap. 
This  is  the  least  amount  of  overlap  allowed  for  the  GLOT. 

An  example  of  the  spectral  coefficients  calculated  in  the  forward  transform  is  given  in 
Figure  4.8.  The  coefficients  in  this  interval  correspond  to  input  samples  from  4096  to  4608. 
Figure  4.7  demonstrates  that  this  interval  contains  significant  speech  data,  yet  our  spectral 
coefficients  seem  to  be  insignificant  in  the  high  frequency  range.  Although  the  exact  nature 
of  the  response  of  the  system  to  eliminating  coefficients  is  unknown  until  tested,  this  initial 
plot  lends  confidence  to  the  idea  that  a  significant  amount  of  compression  may  be  possible 
using  the  GLOT,  without  significant  degradation  of  the  signal  and  without  extremely 
burdensome  computational  complexity.  First,  however,  it  was  necessary  to  determine  the 
ability  of  the  algorithm  to  reconstruct  the  signal  using  all  spectral  components. 

The  signal  was  reconstructed  using  all  of  the  spectral  coefficients  and  the  result  is 
given  in  Figure  4.9.  If  this  plot  looks  quite  similar  to  that  in  Figure  4.7,  it  should.  They  are 
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Figure  4.8.  Sample  Spectral  Coefficients  in  One  Window  on  Interval  [4096,4608] 

exactly  the  same.  The  ESPS  program  used  to  generate  recordings  required  integer  inputs; 
the  output  signal,  therefore,  was  rounded  to  the  nearest  integer  before  being  sent  to  a  file. 
The  input  signal  also  consisted  of  integers.  After  reconstructing  this  speech  signal  using 
all  the  coefficients  and  a  one  point  overlap,  each  of  the  60,080  samples  of  the  output  and 
input  were  exactly  equivalent. 

The  perfect  reconstruction  of  the  first  and  last  data  segments,  using  the  method 
outlined  in  Section  2.2,  was  also  confirmed.  Using  the  GLOT  in  this  manner  provides 
for  a  solution  to  finite  length  signal  processing  previously  unavailable.  This  problem  is  of 
significance  in  applications  such  as  image  processing,  and  of  lesser  importance  in  speech 
processing  where  data  flow  is  continuous.  This  result  demonstrates  the  potential  appli¬ 
cation  of  the  GLOT  to  a  wide  range  of  signal  processing  problems  beyond  the  speech 
application  discussed  later  in  this  section.  The  three  goals  of  a  nonexpansive  transform, 
perfect  reconstruction,  and  no  additional  complexity  were  all  met  using  this  approach. 

As  a  separate  effort,  the  independent  interval  processing  program  detailed  in  Ap¬ 
pendix  C  was  tested.  This  program  provided  perfect  reconstruction  of  the  input  signal 
using  512  sample  intervals  and  a  10%  overlap.  As  this  concept  was  developed  late  in  the 
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Figure  4.9.  Reconstructed  Speech  Signal 

thesis  effort,  no  further  experiments  were  performed  using  this  program.  Chapter  5  dis¬ 
cusses  the  potential  for  follow-on  research  in  the  area  of  independent  interval  processing 
using  the  GLOT. 

Obtaining  perfect  reconstruction  with  a  one  point  overlap  was  a  very  significant 
result.  The  previous  work  in  this  area,  in  particular  by  Malvar  (21),  (18)  and  Princen 
and  Bradley  (27),  were  based  on  a  50%  overlap  which  requires  more  processing  since  each 
point  is  processed  as  part  of  two  intervals.  Further,  Malvar ’s  fast  algorithm  approximates 
the  LOT  (21:556-557),  while  the  GLOT  requires  no  approximations.  Since  the  GLOT  is 
of  the  same  complexity  as  the  FFT,  it  is  comparable  to  Malvar’s  algorithm  (21:553).  This 
result  also  formed  a  baseline  from  which  compression  techniques  could  be  investigated  in 
an  attempt  to  acquire  the  basic  elements  of  a  low  bit  rate  speech  encoder.  As  discussed  in 
Section  1.2,  perfect  reconstruction,  when  using  all  coefficients,  is  an  essential  requirement  to 
further  pursuing  the  merits  of  a  specific  decomposition/reconstruction  algorithm.  Because 
of  the  excellent  results  with  a  one  point  overlap,  subsequent  efforts  focussing  on  compression 
techniques  were  primarily  based  on  this  minimal  overlap. 
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4-6  Exploring  Various  Compression  Techniques 

Building  on  the  success  of  the  previous  section,  various  compression  techniques  could 
be  investigated  in  an  effort  to  reduce  the  number  of  bits  required  to  represent  a  speech 
signal  with  a  high  degree  of  accuracy.  One  problem  associated  with  determining  the  most 
promising  approaches  to  compressing  speech  is  that  success  is  defined  by  the  ear.  Firm 
quantitative  measures  for  speech  coding  success  have  yet  to  be  developed,  and  subjective 
measures  are  often  employed  to  quantify  a  degree  of  success.  Ongoing  research  attempts 
to  link  objective  measures  with  the  subjective  success  or  failure  of  various  compression 
techniques  (see,  for  example,  Quackenbush  et  al.  (28)).  No  firm  guidelines  could  be  used 
to  determine  which  approaches  would  be  most  promising,  because  the  approach  best  suited 
to  any  algorithm  when  dealing  with  speech  signals  is  idiosyncratic  to  the  method  employed 
(15).  The  following  sections  outline  the  compression  techniques  investigated  pursuant  to 
this  thesis.  Each  of  these  techniques  required  some  processing  in  addition  to  that  accounted 
for  in  Appendix  A.  Additional  code  used  in  conjunction  with  the  following  techniques  will 
be  provided  as  part  of  the  appropriate  section. 

4.6.1  Absolute  Spectral  Thresholding  Absolute  Thresholding  was  one  of  the  most 
basic  and  effective  techniques  employed  in  reducing  the  number  of  spectral  coefficients, 
while  maintaining  a  high  degree  of  speech  quality.  The  code  given  in  Figure  4.10  was 
added  to  Procedure  Do  Work  following  the  procedure  call  to  Compute  Coefficients  and 
replaces  the  print  loop  at  that  location. 

As  noted  in  Section  1.4,  the  three  goals  of  a  high  quality,  high  compression,  low 
complexity  coder  must  all  be  weighed.  This  technique  has  a  relatively  low  complexity, 
since  only  a  simple  threshold  comparison  is  required  beyond  the  basic  transform.  While 
setting  a  threshold  is  a  simple  technique,  the  choice  of  a  proper  threshold  is  arbitrary. 
In  general,  the  proper  threshold  will  depend  on  the  data  being  processed  and  effect  the 
compression  capability  of  the  system.  With  this  in  mind,  testing  was  accomplished  to 
determine  the  quality  versus  compression  trade-off.  Four  different  compression  ratios  were 
used  to  process  the  test  phrase,  corresponding  to  85.0%,  92.5%, 96.2%,  and  98.1%.  This 
means  that  the  stated  percentage  of  least  significant  coefficients  were  not  used,  or  “zeroed 
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for  X  in  St  art. Dat  a  . ( End-Data- 1)  loop 

if  abs(Transfornied-Data(x))  <  Alpha-Threshold  then 
Transfornied-Data(x)  :=  0.0; 
if  j  <  (Partitions'last  -  1)  then 

Alpha.Compression-Counter  :=  Alpha-Coinpressioii-Counter  +  1; 
end  if; 
end  if; 

Float  Jo. put(Alpharile, Transformed  .Data(x)); 

Text  Jo.newJine(  Alphafile); 
end  loop; 


F  igure  4.10.  Additional  Code  for  the  Absolute  Threshold  Compression  Technique 


out”  prior  to  reconstruction  of  the  signal.  The  following  assumptions  allowed  for  a  bit  rate 
calculation  corresponding  to  the  percentages  chosen  above  : 

1.  Downsampling  the  input  from  16  kHz  to  8  kHz  would  not  significantly  affect  the 
results. 

2.  The  spectral  coefficients  could  be  represented  with  4  bits. 

3.  Additional  overhead  bits  need  not  be  accounted  for. 


Bit  rates  were  then  computed  by  : 


„,bits,  _ data.sample,  coefficient  .  bits  , 

Rhi - )  =  8000( - )  ♦  1(-T-- - j-)  ♦  4( — — — -)  ♦  (1  -  C).  (4.2) 

sec  sec  data  sample  cocthcient 


where  is  the  bit  rate  of  the  system,  and  C  is  the  compression  ratio  of  eliminated 
coefficients.  This  formula  yields  bit  rates  of  4800,2400, 1200,  and  600  bits/sec  when  using 
the  compression  ratios  mentioned  above.  These  bit  rates  are  commonly  used  in  speech 
processing  systems.  The  same  data  rates  were  tested  here  to  provide  a  bit  rate  comparison 
with  existing  techniques. 

The  thresholded  coefficients  of  one  interval  for  compression  percentages  of  85.0%  and 
98.1%  are  shown  in  Figure  4.11  and  Figure  4.12,  respectively.  Notice  the  significant  re¬ 
duction  in  the  number  of  non-zeroed  coefficients,  especially  in  Figure  4.12,  when  compared 
to  the  unthresholded  coefficients  in  Figure  4.8.  The  first  10,000  points  of  the  output  wave- 
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Figure  4.11.  Absolute  Threshold  Spectral  Coefficients  (85.0%)  on  Interval  [4096,4608] 
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Figure  4.12.  Absolute  Threshold  Spectral  Coefficients  (98.1%)  on  Interval  [4096,4608] 
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Figure  4.13.  Absolute  Threshold  Output  Waveform  with  85.0%  Compression 


form  for  these  same  compression  percentages  are  given  in  Figures  4.13  and  4.14.  When 
compared  to  Figure  4.9,  Figure  4.13  provides  a  better  match  to  the  perfectly  reconstructed 
waveform  than  Figure  4.14,  especially  when  the  speech  is  of  low  amplitude.  Of  course, 
the  final  test  of  the  success  or  failure  of  these  outputs  in  matching  the  input  waveform 
depends  on  the  ear.  After  listening  to  the  absolute  thresholded  outputs,  some  interesting 
observations  were  made.  The  reconstructed  signal  at  an  85.0%  threshold  level  was  of 
high  quality.  In  each  recording,  noise  was  evident  in  the  form  of  “pinging”  or  extraneous 
tones.  The  speaker’s  voice,  however,  especially  at  85.0%  thresholding,  was  not  distorted. 
In  other  words,  the  speaker,  as  well  as  the  speech,  could  be  discerned  in  the  recording.  In 
addition  to  the  background  tones,  the  other  form  of  noise  was  the  elimination  of  parts  of 
speech  near  word  boundaries.  This  can  be  seen  in  Figure  4.14  as  the  degraded  reconstruc¬ 
tion  just  after  sample  2000  and  again  between  sample  6000  and  8000.  Even  at  the  85.0% 
level,  the  final  syllable  of  the  final  word  “nihilistic”  was  not  discernible,  but  the  problem 
was  most  significant  at  higher  compression  levels.  The  98.1%  thresholded  output  was  not 
recognizable,  due  to  both  forms  of  noise  discussed  above. 
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Figure  4.14.  Absolute  Threshold  Output  Waveform  with  98.1%  Compression 

A  modification  to  this  approach  was  also  implemented  after  viewing  the  results  of  the 
absolute  threshold  technique.  At  high  compression  ratios,  a  single  (or  small  number  of) 
remaining  coefficient(s)  might  result  in  a  reconstruction  containing  sinusoidal  noise  in  the 
entire  interval.  Figure  4.15  is  an  example  of  two  such  intervals.  In  the  two  adjacent  intervals 
shown,  only  one  coefficient  was  of  large  enough  magnitude  to  survive  the  thresholding. 
The  resulting  reconstruction  is  a  low  frequency  sinusoid  from  the  low  frequency  coefficient 
in  the  left  interval,  and  a  high  frequency  sinusoid  from  the  high  frequency  coefficient  in 
the  right  interval.  Figure  4.16  yields  this  reconstruction,  which  is  a  form  of  noise  in 
the  output  wavefoii  i.  (Note  also  the  180 deg  phase  shift  in  the  low  frequency  sinusoid 
resulting  from  a  negative  coefficient,  as  expected).  To  lessen  this  problem,  and  allow 
a  more  judicious  choice  of  important  coefficients,  the  following  rule  was  applied  :  any 
coefficient  above  the  threshold,  but  with  neither  adjacent  coefficient  above  the  threshold, 
will  be  zeroed  as  if  it  had  not  met  the  threshold  criterion.  Besides  the  visible  sinusoidal 
noise,  an  underlying  a.ssumption  in  making  this  rule  was  that  coefficients  containing  a 
significant  amount  of  signal  information  are  generally  grouped  together.  This  assumption 
was  based  on  viewing  numerous  spectral  plots.  Based  on  this  assumption,  the  converse  of 
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Figure  4.15.  Single  Spectral  Coefficient  in  Each  of  Two  Adjacent  Intervals 
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Figure  4.16.  Reconstructed  Data  in  Two  Adjacent  Windows,  [38400—  39424] 
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the  aforementioned  rule  was  also  applied  :  any  coefficient  below  the  threshold,  but  with 
both  adjacent  coefficients  above  the  threshold,  will  be  retained  for  reconstruction.  The 
code  for  this  modified  absolute  threshold  method  is  given  in  Figure  4.17.  Recordings  of  the 
output  waveforms  using  this  modified  absolute  threshold  approach  were  indistinguishable 
from  the  recordings  when  using  the  original  approach.  The  comparison  was  performed  for 
each  of  the  four  compression  percentages.  Some  differences  in  the  waveform  were  apparent, 
but  they  were  not  significant  enough  to  be  distinguished  by  the  ear. 

4.6.2  Relative  Spectral  Thresholding  The  next  approach  was  similar  to  the  absolute 
threshold  method,  but  each  coefficient  was  scaled  prior  to  comparison  with  the  threshold. 
The  scale  factor  was  simply  the  sum  of  all  coefficients  in  the  window,  or  the  integral 
of  coefficients  in  the  window.  Adding  this  parameter  into  the  calculations  was  an  at¬ 
tempt  at  providing  a  better  reconstruction  for  intervals  containing  edges  of  words  or  near 
speech/non-speech  boundaries.  This  portion  of  the  reconstructed  signal  was  degraded 
substantially  using  the  absolute  threshold  method,  especially  at  higher  comp’-ession  per¬ 
centages  (or  lower  bit  rates),  as  discussed  in  Section  4.6.1.  The  code  for  this  new  approach 
is  given  in  Figure  4.18.  Figures  4.19  and  4.20  are  the  reconstructed  waveform  after  ap¬ 
plying  the  relative  threshold  technique.  The  85.0%  reconstruction  has  some  spikes  in  the 
output  at  various  sample  locations.  At  the  93.1%  level,  the  output  is  severely  distorted  in 
many  intervals.  When  comparing  a  recording  of  these  outputs  to  the  absolute  threshold 
technique  for  similar  compression  percentages,  the  absolute  threshold  method  produced  a 
better  quality  reproduction.  There  was  more  background  noise  or  “tones"’  in  the  sound 
of  the  relative  thresholded  reproduction,  and  the  sound  of  the  voice  also  changed  slightly 
(affecting  the  recognition).  The  93.1%  compressed  reconstruction  was  mostly  pure  noise 
to  the  ear,  as  would  be  expected  from  viewing  Figure  4.20.  At  85.0%  compression,  the 
relative  threshold  signal,  while  worse  than  its  absolute  threshold  equivalent  in  general, 
had  one  redeeming  quality.  The  signal  at  word  boundaries  was  often  reproduced  using 
the  relative  threshold  technique  whereas  it  was  eliminated  in  the  absolute  threshold  out¬ 
put.  The  relative  threshold  method  reproduced  the  “...tic”  of  “nihilistic”,  discussed  in 
Section  4.6.1.  This  advan.  ige  is  the  result  of  forcing  coefficients  from  each  window  to  be 
retained,  thereby  allowing  a  reproduction  even  where  there  is  minimal  signal  energy.  The 
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— Take  care  of  first  coefficient  separately 
— so  the  first  time  through  the  following 
— loop  there  is  a  (x  —  1)  point  to  look  at 
if  abs(Transformed_Data(Start_Data))  <  Alpha_Threshold  then 
Transformed_Data(Start-Data)  :=  0.0; 
if  j  <  (Partitions’last -1)  then 

Alpha.Compression.Counter  :=  Alpha.Compression.Counter  +  1; 
end  if; 
end  if; 

Float  Jo. put(Alphafile, Transformed  JData(Start_Data)); 
TextJo.New.Line(Alphafile); 

for  X  in  (Start.Data+l)..(End-Data-l)  loop 

if  abs(Transformed-Data(x))  <  Alpha.Threshold  then 
if  Transformed_Data(x-l)  =  0.0  then 
TransformedJDatatx)  :=  0.0; 
ifj<  (Partitions’last -1)  then 

Alpha.Compression.Counter  :=  Alpha.Compression.Counter  +  1; 
end  if; 

elsif  abs(Transformed.Data(x+l))  <  Alpha.Threshold  then 
Transformed  J)ata(x)  :=  0.0; 
ifj  <  (Partitions’last -1)  then 

Alpha.Compression.Counter  :=  Alpha.Compression.Counter  +  1; 
end  if; 
end  if; 

elsif  Transformed.Data(x-l)  =  0.0  then 

if  abs(Transformed.Data(x+l))  <  Alpha.Threshold  then 
Transformed J)ata(x)  :=  0.0; 
if  j  <  (Partitions’last -1)  then 

Alpha.Compression.Counter  Alpha.Compression.Counter  +  1; 
end  if; 
end  if; 
end  if; 

Float  Jo.  put(Alphafile,'IVansformedJ)ata(x)); 

TextJo.New_Line(Alphafile); 
end  loop; 


Figure  4.17. 


Additional  Code  for  the  Modified  Absolute  Threshold  Compression 
Technique 
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Alpha-Sum  ;=  0.0; 

—  Alpha-Sum  is  recomputed  for  each  interval 

for  X  in  Start-Data. .(End-Data-1)  loop 

Alpha-Sum  :=  Transformed-Data(x)  +  Alpha-Sum; 
end  loop; 

for  X  in  Start-Data.. (End-Data- 1)  loop 

if  (Transformed-Data(x)/Alpha_Sum)  <  ReLThreshold  then 
Transformed-Data(x)  :=  0.0; 
if  j<  (Partitions’last  -1)  then 

Alpha-Compression-Counter  ;=  Alpha-Compression-Counter  -f  1; 
end  if; 
end  if; 

FloatJo.Put(Alphafile,Transformed-Data(x)); 
Text-Io.NewXine(Alphafile); 
end  loop; 

Figure  4.18.  Additional  Code  for  the  Relative  Thresholding  Compression  Technique 
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Figure  4.19.  Relative  Threshold  Output  Waveform  with  85.0%  Compression 
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Figure  4.20.  Relative  Threshold  Output  Waveform  with  98.1%  Compression 

tradeoff  was  that  coefficients  were  removed  from  the  windows  containing  a  strong  signal, 
creating  excess  noise  in  numerous  intervals.  This  noise  overshadowed  any  benefit  resulting 
from  use  of  the  relative  threshold  technique. 

4-6.3  Polynomial  Fitting  of  the  Spectrum  All  of  the  speech  compression  techniques 
discussed  thus  far  had  unique  merits  and  shortcomings,  but  none  was  far  superior  in  per¬ 
formance  to  the  others.  The  simplest  technique,  the  basic  absolute  threshold,  provided  the 
best  reconstruction,  yet  none  of  the  others  were  significantly  more  complex.  The  amount 
of  required  coefficient  processing  was  an  important  parameter  to  minimize.  For  this  rea¬ 
son,  the  initial  methods  tested  were  the  least  complex.  It  became  apparent,  however,  that 
further  reductions  in  the  number  of  coefficients  required  for  high  quality  speech  recon¬ 
struction  would  require  more  complicated  spectral  processing.  Viewing  the  unthresholded 
coefficients,  it  appeared  is  if  they  might  be  modelled  by  a  polynomial  curve.  The  follow¬ 
ing  rules  were  established  as  a  criteria  for  what  information  would  be  extracted  from  the 
coefficients  and  used  in  the  reconstruction  : 

1.  A  noise  threshold  will  be  established  below  which  no  coefficients  will  be  kept. 
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Figure  4.21.  Sample  Spectral  Coefficients  Tagged  for  Polynomial  Fit  Processing 

2.  From  the  remaining  coefficients,  the  local  maxima  and  minima  for  both  the  positive 
and  negative  coefficients  will  be  retained. 

3.  Of  those  coefficients  above  the  noise  threshold,  but  not  maintained  as  a  maxima  or 
minima,  the  phase  will  be  maintained. 

Figure  4.21  is  a  portion  of  the  coefficients  given  by  Figure  4.8.  The  first  three  mtixima 
or  minima  which  would  be  saved  in  this  interval  based  on  the  above  rules  are  labeled 
appropriately.  The  noise  threshold  was  assumed  to  be  below  any  of  the  coefficients  in 
this  figure.  The  following  rules  governed  reconstruction  of  the  signal  given  this  saved 
information  : 

1.  A  deterministic  polynomial  curve  will  be  fitted  using  the  strictly  positive  or  strictly 
negative  local  maxima,  local  minima,  and  local  maxima  (see  Figure  4.21)  as  the  left, 
middle,  and  right  points  on  the  curve  respectively.  This  process  will  result  in  both 
a  positive  and  negative  polynomial  curve  being  defined  at  any  given  point  in  the 
spectrum. 
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2.  Any  coefficient  which  was  not  a  maxima  or  minima,  but  was  above  the  noise  thresh¬ 
old,  should  be  estimated  by  the  value  of  the  appropriately  phased  polynomial  curve 
defined  at  that  location. 

3.  Any  coefficient  which  was  below  the  noise  threshold  will  be  assumed  to  be  zero. 

This  method  required  the  addition  of  three  new  procedures  outlined  in  Appendix  B. 
Procedure  Maxima-Minima  extracted  the  positive  and  negative  maxima  and  minima,  the 
required  phase  information,  and  some  initial  conditions  for  use  in  the  reconstruction.  Pro¬ 
cedure  PolynomiaLFit  defined  the  parameters  to  be  used  in  procedure  FiLSlice,  which 
computed  the  actual  polynomial  and  reconstructs  the  coefficients.  At  this  point,  the  basic 
reconstruction  is  performed  on  the  estimated  coefficients. 

This  method  is  different  from  any  of  the  other  methods  in  that  the  amount  of  pos¬ 
sible  compression  is  set  by  the  data  rather  than  varied  based  on  a  compression  threshold. 
Although  some  range  of  compression  can  be  achieved  by  modifying  the  noise  threshold, 
the  compression  percentage  is  essentially  the  ratio  of  relative  maxima  and  minima  spectral 
coefficients  to  data  samples.  An  advantage  of  this  technique  is  that  the  choice  of  threshold 
is  not  arbitrary,  but  related  only  to  the  amount  of  noise  in  the  signal.  Viewing  the  data 
in  non-speech  intervals,  it  was  determined  that  the  noise  reached  an  amplitude  level  of  ap¬ 
proximately  35.  Setting  the  threshold  to  this  level  resulted  in  a  compression  ratio  of  20%. 
A  threshold  of  60  yielded  a  compression  ratio  of  85.0%,  analogous  to  data  accumulated  for 
previous  methods.  Varying  the  threshold  from  35  to  60  did  not  perceptibly  alter  either  the 
viewed  or  recorded  reconstructed  output.  The  coefficients  remaining  after  extracting  the 
maxima  and  minima  for  one  interval  are  given  in  Figure  4.22.  The  coefficients  in  this  inter¬ 
val  are  much  more  sparse  than  the  original  unthresholded  coefficients  given  in  Figure  4.8. 
The  reconstructed  output  with  85.0%  compression  using  the  polynomial  fit  technique  is 
given  in  Figure  4.23. 

Figures  4.24  through  4.26  are  provided  to  demonstrate  the  most  important  aspect  of 
the  polynomial  fit  technique,  and  the  step  which  creates  noise  in  the  reconstructed  output. 
Figure  4.24  is  a  plot  of  unthresholded  spectral  coefficients  64  -  128  from  the  interval 
analyzed  in  Figure  4.8.  Figure  4.25  are  the  remaining  (positive  and  negative)  maxima  and 
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Figure  4.22.  Polynomial  Fit  Spectral  Coefficients  (85.0%)  on  Interval  [4096,4608] 


Figure  4.23.  Polynomial  Fit  Output  Waveform  with  85.0%  Compression 
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Figure  4.24.  Unthresholded  Coefficients  64  -  128  of  Figure  4.8 


128 


minima  which  will  be  used  to  reconstruct  the  signal,  along  with  a  sign  vector.  Figure  4.26 
gives  the  reconstructed  or  estimated  spectral  coefficients  based  on  the  polynomial  fit  criteria 
discussed  above.  Comparing  Figures  4.24  and  4.26,  it  is  apparent  that  in  some  regions, 
the  polynomial  fit  accurately  estimates  the  actual  coefficients.  In  other  regions,  however, 
the  coefficients  cannot  be  modelled  accurately  with  a  polynomial  curve.  Reconstructing 
the  signal  using  coefficients  which  have  not  been  accurately  modelled  creates  noise  in 
the  output  waveform.  The  full  extent  of  this  noise  can  only  be  ascertained,  however,  by 
listening  to  the  recorded  output. 

The  polynomial  fit  output  was  of  the  same  audio  quality  as  the  absolute  threshold 
output  of  similar  compression  ratio.  It  also  had  the  advantage,  as  in  the  case  of  the  relative 
threshold  technique,  of  retaining  some  speech  sounds  at  word  boundaries  which  were  elim¬ 
inated  in  the  absolute  threshold  output.  In  this  sense,  the  polynomial  fit  technique  was 
marginally  superior  in  audio  reproduction,  at  the  expense  of  some  additional  computations. 
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Figure  4.25.  Polynomial  Fit  Spectral  Coefficients  64  -  128  of  Figure  4.22 
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Figure  4.26.  Reconstructed  Polynomial  Fit  Spectral  Coefficients  64  -  128  of  Figure  4.25 
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4.6.4  Cepstral  Processing  Homomorphic  filtering,  or  homomorphic  deconvolution, 
is  a  process  designed  to  separate  convolved  components  of  a  signal  (25:807).  Cepstral  pro¬ 
cessing  is  a  specific  homomorphic  filtering  technique  with  numerous  applications,  including 
speech  (6:1428).  The  term  “cepstrum”  comes  from  the  interchange  of  each  letter  in  the 
first  syllable  of  “spectrum”.  Using  this  same  naming  convention,  terms  such  as  quefrency 
and  liftering  correspond  to  frequency  and  filtering  in  the  cepstral  domain.  The  power  cep¬ 
strum  is  generally  defined  as  the  power  spectrum  of  the  logarithm  of  the  power  spectrum 
of  the  function,  and  is  used  in  detection  or  estimation  of  signal  parameters  (6:1429).  An¬ 
other  technique,  the  complex  cepstrum,  retains  the  phase  information  of  the  data,  and  can 
therefore  be  used  in  reconstruction  applications  (6:1430).  The  equation  for  the  complex 
cepstrum  (25:782)  is  given  in  Equation  4.3  below: 

Cm  =  ^  y  log  I  5(e’“')  I  (4.3) 

The  disadvantage  of  this  technique  is  the  complexity  involved  in  calculating  the  complex 
logarithm.  The  transform  used  in  this  type  of  system  is  a  Fourier  transform.  The  technique 
discussed  in  this  section  involves  application  of  cepstral-like  techniques  using  the  GLOT. 

The  coefficients  of  the  GLOT  are  real  valued,  yet  in  magnitude  and  phase  represen¬ 
tation,  they  contain  a  phase  term  of  either  zero  or  180  deg  which  must  be  maintained  to 
reconstruct  the  signal.  Recognizing  that  the  phase  term  needed  to  be  maintained,  but  that 
it  could  only  be  one  of  two  specific  values,  the  following  definition  was  derived  to  govern  a 
cepstral-like  transformation  as  it  relates  to  the  GLOT  : 


Cj.k(x)  =  ^F[ln  I  (£?p(.s(x_,,t))  Ij.  (4.4) 

where  s(xj  *)  is  the  data,  Qp  is  the  forward  GLOT,  and  Cj  k{x)  are  the  resulting  cepstral-like 
coefficients.  This  new  transform  actually  encompasses  parts  of  both  the  power  cepstrum 
and  complex  cepstrum  described  earlier.  As  in  the  power  cepstrum,  a  complex  logarithm  is 
not  required.  The  magnitude  operation  in  equation  4.4  facilitates  this  simplified  operation. 
Because  the  phase  information  can  only  be  one  of  two  values,  it  can  easily  be  extracted  as 
part  of  the  algorithm  and  used  in  the  reconstruction.  The  phase  information  is  preserved 
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for  reconstruction,  as  in  the  complex  cepstrum  given  in  Equation  4.3,  without  using  the 
complex  logarithm.  By  tailoring  cepstral-type  operations  specifically  to  the  GLOT,  the 
best  of  both  techniques  can  be  exploited.  Given  cepstral-like  coefficients  defined  by  £j  k{x), 
the  signal  can  be  reconstructed  by  : 

Sj.k  =  Gi[Pi.kGi(CjA^))]-  (4.5) 

where  Qi  is  the  inverse  GLOT,  Pj  ^  is  a  phase  vector  correcting  for  the  magnitude  scaling 
in  equation  4.4,  and  5j  *  is  the  reconstructed  signal  such  that  *  «  s{xj  ie). 

Using  this  new  decomposition  and  reconstruction  of  the  signal,  various  compression 
techniques  of  the  cepstral  coefficients  (mocpression  techniques  ?)  can  be  tested.  If  this 
new  transform  can  “pack”  more  signal  information  into  fewer  coefficients,  then  greater 
compression  can  result  by  retaining  only  the  information  packed  cepstral  coefficients. 

This  transformation  is  different  from  the  original  GLOT.  Therefore,  the  first  test  was 
to  use  all  cepstral  coefficients  in  reconstruction  of  the  signal.  The  cepstral  coefficients  in 
one  window  are  given  in  Figure  4.27.  The  output  signal  was  an  exact  copy  of  the  input 
signal  under  these  conditions.  In  this  test,  the  perfect  reconstruction  criteria  was  satisfied, 
and  correct  implementation  of  the  theory  outlined  above  was  also  validated.  The  next  two 
sections  detail  two  methods  employed  in  an  effort  to  determine  a  high  quality,  low  bit  rate 
technique  for  speech  processing. 

4. 6.4-1  Cepstral  Liftering  The  cepstral  coefficients  near  the  origin  (low  que- 
frency)  define  the  spectral  envelope  (26:203-204).  Based  on  this  observation,  one  technique 
which  can  be  employed  is  to  retain  the  first  x  coefficients  to  be  used  in  the  reconstruction 
of  the  signal,  where  x  is  can  vary  from  one  to  the  total  number  of  cepstral  coefficients. 
The  compression  ratio  is  then  given  by  x/iV,  where  N  is  the  number  of  samples  (and 
cepstral  coefficients)  in  each  interval.  This  technique  involves  liftering  (filtering)  the  low 
quefrency  cepstral  coefficients.  Figure  4.27  displays  low  quefrency  coefficients  having  higher 
amplitude  than  the  rest  of  the  cepstral  coefficients,  and  liftering  will  pass  these  important 
coefficients.  The  same  four  compression  percentages  used  earlier  were  tested  here.  The  out- 
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Figure  4.27.  Non-liftered  Non-thresholded  Cepstral  Coefficients  in  One  Window  on  In¬ 
terval  [4096,4608] 

puts  for  85.0%  and  98.1%  compression  of  the  cepstral  coefficients  are  given  in  Figures  4.28 
and  4.29,  respectively. 

In  all  previous  plots,  higher  quality  recorded  outputs  had  waveforms  that  matched  the 
input  waveform  to  a  high  degree.  Conversely,  output  waveforms  which  did  not  match  the 
input  waveform  were  of  poor  audio  quality.  This  trend  does  not  indicate,  however,  that  a 
correlation  always  exists  between  reconstructed  speech  quality  and  waveshape.  Comparing 
Figures  4.28  and  4.29  indicates  strong  similarity  between  the  outputs.  When  listening  to 
these  two  outputs,  however,  the  output  in  Figure  4.28  was  of  significantly  higher  audio 
quality. 

The  cepstral  lifter  technique  distorted  the  speakers  voice  slightly  at  the  85.0%  com¬ 
pression  level.  There  were  no  extraneous  tones  in  the  output  as  in  the  case  of  the  absolute 
threshold  output.  The  cepstral  liftering  seemed  to  garble  the  speaker’s  voice,  and  it  was 
not  as  clear  as  the  absolute  threshold  output  at  this  compression  level.  At  higher  com¬ 
pression  ratios,  the  cepstral  lifter  technique  provided  much  higher  quality  outputs  when 
compared  to  the  absolute  threshold  method.  The  transition  from  85.0%  to  92.5%  did  not 
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Figure  4.28.  Output  Waveform  After  Cepstral  Liftering  and  Performing  85.0% 
Compression 


Figure  4.29.  Output  Waveform  After  Cepstral  Liftering  and  Performing  98.1% 
Compression 
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Figure  4.30.  Output  Waveform  After  Cepstral  Liftering  and  Performing  99.1% 
Compression 

audibly  detract  from  the  quality  of  the  reconstructed  output,  and  further  compression  to 
96.2%  and  98.1%  were  the  highest  quality  reproductions  at  these  levels  compared  to  any 
of  the  previous  techniques.  This  success  prompted  another  transformation  for  a  bit  rate 
of  300  bit/set  as  computed  in  Section  4.6.1.  300  bits/sec  corresponds  to  using  less  than 
one  percent  of  the  coefR<  ients  in  reconstruction  of  the  signal!  Figure  4.30  is  a  plot  of  the 
output  resulting  from  this  reconstruction.  The  speech  was  recognizable  at  this  level,  but 
contained  significant  noise.  The  waveform  is  very  similar  to  the  waveforms  in  Figures  4.28 
and  4.29,  although  the  audio  quality  of  all  three  is  different. 

This  technique  was  different  from  all  other  techniques  in  that  the  number  of  coef¬ 
ficients  retained  in  each  window  remained  constant.  Less  overhead  would  be  required  to 
pass  the  necessary  reconstruction  information  through  a  communications  channel,  since  the 
location  and  number  of  coefficients  in  each  window  is  preordained.  In  all  other  schemes, 
the  number  of  coefficients  and  their  unique  position  in  frequency  would  need  to  be  passed 
in  addition  to  the  value  of  the  coefficients.  The  disadvantage  of  this  technique  is  that  the 
same  number  of  coefficients  are  used  to  represent  windows  containing  no  speech  as  those 
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Figure  4.31.  Thresholdcd  Cepstral  Coefficients  (85.0%)  in  One  Window  on  Interval 
[4096,4608] 


containing  speech  in  the  entire  interval.  If  12  coefficients  are  liftered  from  each  window, 
and  there  are  10  windows  where  no  speech  is  present,  then  120  cepstral  coefficients  are 
retained  which  have  no  information.  To  overcome  this  indiscriminate  coefficient  retention, 
the  technique  of  the  next  section  was  investigated. 

4- 6. 4 -2  Cepstral  Thresholding  Rather  than  liftering  the  low  quefrency  cep¬ 
stral  coefficients  from  each  window,  a  threshold  was  used  to  retain  a  percentage  of  the 
highest  amplitude  coefficients.  Because  the  low  quefrency  coefficients  are  genera'Iy  of 
higher  amplitude  (see  Figure  4.27),  this  technique  has  the  effect  of  low  pass  liftering.  The 
difference  here  is  that  more  low  quefrency  cepstral  coefficients  will  be  retained  in  speech 
intervals  and  less  will  be  retained  in  non  speech  intervals.  This  technique  is  completely 
analogous  to  the  spectral  thresholding  technique  of  Section  4.6.1.  The  coefficients  re¬ 
tained  i!i  one  example  window  for  compression  percentages  85.0%  and  98.1%  are  given 
in  Figures  4.31  and  4.32,  respectively.  Notice  that  some  higher  quefrency  coefficients 
are  retained  in  Figure  4.31,  due  to  the  higher  cepstral  threshold.  When  the  threshold  is 
increased  to  achieve  98.1%  compression,  however,  the  result  is  a  basically  a  low  pass  lifter 
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Figure  4.32.  Thresholded  Cepstral  Coefficients  (98.1%)  in  One  Window  on  Interval 
[4096,4608] 
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of  the  first  12  cepstral  coefficients.  Figure  4.33  gives  only  the  low  quefrency  coefficients  of 
Figure  4.33.  The  only  difference  between  a  12  coefficient  low  pass  lifter  and  the  technique 
resulting  in  Figure  4.33  is  that  one  of  the  12  low  quefrency  coefficients  did  not  meet  the 
threshold,  and  was  zeroed  out  along  with  the  higher  quefrency  coefficients.  The  output 
resulting  from  the  85.0%  compression  level,  given  in  Figures  4.34,  was  closer  in  form  to  the 
input  than  any  of  the  liftering  outputs.  At  higher  compression  ratios,  however,  the  output 
waveforms  were  very  similar  to  the  liftering  output  waveforms  given  above. 

This  technique  detracted  from  the  speaker’s  voice  quality  less  than  the  liftering  tech¬ 
nique.  At  the  85.0%  compression  level,  however,  the  absolute  threshold  technique  still 
had  the  best  quality  reproduction.  At  higher  compression  ratios  (above  98%),  cepstral 
thresholding  outperformed  all  other  techniques  tested. 

Because  of  the  success  of  this  technique  at  high  compression  ratios,  the  99.1%  com¬ 
pression  level  was  tested  as  in  Section  4.6.4. 1  above.  Figures  4.35  and  4.36  give  the  thresh¬ 
olded  cepstral  coefficients  of  one  window  and  resulting  output  waveform  for  this  compres¬ 
sion  level.  Only  the  first  64  cepstral  coefficients  are  shown  in  Figure  4.35,  since  all  others 
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Figure  4.33.  First  64  Thresholded  Cepstral  Coefficients  (98.1%)  in  One  Window  on  In¬ 
terval  [4096,4608] 
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Figure  4.34.  Cepstral  Thresholding  Output  Waveform  with  99.1%  Compression 
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Figure  4.35.  First  64  Thresholded  Cepstral  Coefficients  (99.1%)  in  One  Window  on  In 
terval [4096,4608] 


Figure  4.36.  Cepstral  Thresholding  Output  Waveform  with  99.1%  Compression 
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are  zero.  The  audio  quality  at  this  high  compression  level  was  better  than  the  quality  of 
the  analogous  liftering  output.  The  speech  was  understandable,  although  noise  was  present 
which  altered  the  speaker’s  voice. 

4.6.5  Summary  All  of  the  techniques  given  in  Section  4.6  were  analyzed  to  gain  a 
better  understanding  of  the  GLOT,  and  define  potential  areas  which  merit  further  research 
toward  development  of  a  low  bit  rate  speech  encoder.  Due  to  the  amount  of  previous 
research  in  the  area  of  speech  compression,  and  the  complexities  of  speech  processing,  the 
body  of  experiments  performed  in  this  thesis  is  not  sufficient  to  define  such  an  encoder. 
The  following  chapter  will  focus  on  potential  follow-on  topics  toward  this  goal,  as  well  as 
summarize  the  results  of  this  thesis. 
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V.  Conclusion 


5.1  Summary  and  Conclusions 

This  thesis  was  the  incipiency  of  research  regarding  the  generalized  lapped  transform 
and  its  application  to  the  processing  of  speech  signals.  As  such,  the  initial  research  empha¬ 
sis  was  proper  validation  and  implementation  of  the  theory.  Many  areas  hold  promise  for 
future  research,  and  some  of  these  topics  are  outlined  in  the  next  section.  First,  however, 
a  summary  of  the  results  of  this  thesis  will  be  provided  which  have  laid  the  foundation  for 
subsequent  efforts. 

The  first  stage  in  validating  the  theory  behind  the  generalized  lapped  transform 
was  the  development  of  a  prototype  system  using  the  digital  signal  processing  tool,  PC- 
DSP.  Use  of  a  prototype  system,  especially  when  implementing  an  untested  algorithm,  can 
significantly  reduce  the  overall  time  to  develop  a  working  system.  The  use  of  PC-DSP  was 
especially  helpful  in  making  minor  corrections  to  theory,  and  quickly  validating  certain 
theoretical  predictions.  Specifically,  the  use  of  PC-DSP  experimentally  demonstrated  that 
perfect  reconstruction  of  a  signal  was  possible  when  using  all  spectral  coefficients  in  the 
reconstruction  algorithm.  Further,  testing  of  the  eventual  software  program  was  facilitated, 
because  expected  results  from  each  stage  of  the  algorithm  were  known  a  priori. 

The  most  significant  coding  task  and  result  was  the  development  and  testing  of  a 
fast  Fourier  transform  (FFT)  algorithm.  Although  developed  specifically  to  support  the 
generalized  lapped  transform  ,  this  Ada  based  algorithm  is  of  use  in  any  digital  signal 
processing  application.  Most  FFT  routines  are  limited  to  radix  2  application;  however, 
this  general  purpose  algorithm  can  process  data  of  arbitrary  length.  The  timing  analysis 
in  Chapter  3  also  reveals  that  this  Ada  based  routine,  and  the  Ada  language  in  general,  is 
viable  for  digital  signal  processing  applications. 

The  software  program  outlined  in  Appendix  A  provides  a  decomposition/reconstruc¬ 
tion  algorithm  based  on  a  sinusoidal  basis  function  and  windows  with  variable  length  and 
percent  overlap  with  adjacent  windows.  Initial  testing  confirmed  that  perfect  reconstruc¬ 
tion  of  a  signal  was  achievable  when  using  all  spectral  coefficients  in  the  reconstruction, 
thereby  satisfying  the  first  research  goal  outlined  in  Chapter  1.  This  result  assumes  the 
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signal  is  defined  outside  the  interval  of  interest.  Due  to  the  success  of  the  finite  length 
signal  processing  technique  outlined  in  Sections  2.2  and  4. 3. 1.6,  perfect  reconstruction  was 
obtained  for  the  first  and  last  data  segments  as  well,  even  though  the  signal  was  undefined 
on  one  adjacent  interval.  This  technique  is  the  first  nonexpansive  lapped  transform  based 
solution  which  achieves  perfect  reconstruction  of  a  finite  length  signal,  and  has  potential 
Air  Force  applications  in  areas  such  as  image  coding,  where  the  signal  of  interest  is  at  the 
border  of  the  image.  The  second  research  goal  involving  finite  length  signals  was  satisfied 
based  on  this  solution. 

Previous  research  in  the  area  of  lapped  transforms  has  demonstrated  perfect  re¬ 
construction  for  data  defined  outside  the  interval  of  interest  using  a  50/fc%  overlap,  and 
mathematically  derived  conditions  for  a  more  general  size  overlap  (19)  (8).  This  thesis 
implemented,  for  the  first  time,  a  lapped  transform  algorithm  based  on  the  work  of  Suter 
and  Oxley  (33)  which  achieved  perfect  reconstruction  of  a  speech  signal  with  a  one  point 
overlap  between  adjacent  windows.  This  result  is  significant  because  a  transform  based  on 
a  one  point  overlap  requires  half  the  processing  of  a  transform  based  on  a  50%  overlap. 

The  independent  interval  processing  program  (See  Appendix  C)  provides  a  method 
of  achieving  perfect  reconstruction  of  a  signal  without  overlapping  adjacent  intervals.  This 
result  is  an  extension  of  the  finite  length  signal  solution,  and  previous  lapped  transform 
algorithms  have  not  addressed  this  method  of  processing  a  signal.  This  result  developed 
late  in  the  thesis  cycle.  Consequently,  only  minimal  testing  and  experimentation  for  this 
specific  implementation  of  the  generalized  lapped  transform  was  performed.  Despite  lim¬ 
ited  experimentation,  this  processing  method  has  merits  which  warrant  further  research 
(See  Section  5.2.1). 

The  software  program  of  Appendix  A  was  thoroughly  tested  to  ensure  operation 
in  accordance  with  the  theory  outlined  in  Chapter  2.  After  validating  the  theory,  a  set 
of  experiments  were  performed  to  gain  insight  into  the  applicability  of  the  transform  to 
speech  processing,  and  specifically  speech  compression.  These  experiments  were  based  on 
various  spectral  coefficient  processing  techniques,  including  absolute  thresholding,  relative 
thresholding,  polynomial  fitting,  and  cepstral  processing.  The  absolute  thresholding  tech¬ 
nique  involved  setting  a  threshold  to  determine  which  spectral  coefficients  would  be  used 
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in  reconstruction  of  the  signal.  Any  coefficient  of  amplitude  greater  than  the  threshold 
was  retained,  while  any  coefficient  not  meeting  the  threshold  was  set  to  zero.  This  tech¬ 
nique  provided  the  best  results  in  terms  of  audio  quality  at  compression  ratios  up  to  85%, 
although  some  perceptible  noise  was  apparent  in  the  reconstructed  output.  The  relative 
threshold  technique  used  the  same  retention  criteria,  although  the  threshold  was  unique  to 
each  window.  This  technique  had  the  advantage  of  better  reproduction  of  the  beginning  or 
ending  of  words,  but  overall,  was  inferior  to  the  absolute  threshold  technique  in  terms  of 
noise.  The  polynomial  fitting  technique  involved  representation  of  the  spectral  coefficients 
with  a  polynomial  curve  defined  by  the  relative  maxima  and  minima  coefficients.  This 
technique  provided  outputs  equal  in  audio  quality  to  the  absolute  threshold  technique  for 
85%  compression.  In  addition,  the  polynomial  fitting  technique  provided  better  speech 
reproduction  at  word  boundaries  (analogous  to  the  relative  threshold  technique).  This 
technique  was  inherently  limited  to  a  compression  ratio  of  85%.  The  most  promising  tech¬ 
nique  investigated  was  the  cepstral  processing  technique,  which  produced  an  intelligible 
output  when  using  less  than  one  percent  of  the  cepstral  coefficients  in  reconstruction  of 
the  signal.  Two  variations  on  cepstral  processing  were  investigated.  The  first  method  em¬ 
ployed  a  low  pass  lifter  (filter)  to  the  cepstral  coefficients,  retaining  the  first  x  coefficients. 
Varying  x  defined  a  compression  percentage.  The  second  method  was  analogous  to  the 
absolute  threshold  method,  only  in  the  cepstral  domain.  A  percentage  of  cepstral  coeffi¬ 
cients  were  retained  based  on  a  threshold  level.  The  second  cepstral  technique  produced 
a  better  quality  audio  output  compared  to  the  first  method.  The  trade-off  is  that  the 
first  method  requires  overhead  to  send  the  position  of  retained  cepstral  coefficients  in  each 
window.  The  second  method  involves  reconstruction  based  on  a  set  number  and  position 
of  low  quefrency  cepstral  coefficients,  and  only  the  coefficient  amplitude  need  be  encoded. 
.\11  of  the  methods  described  were  investigated  to  determine  the  applicability  of  the  GLOT 
to  speech  processing  as  outlined  in  Goal  three  of  Chapter  1.  Further  research  is  needed 
(see  Section  5. 2. .3)  to  determine  if  these  techniques  can  be  exploited  to  define  a  complete 
low  bit  rate  speech  encoder. 


5.2  Follow-On  Research  Areas 


5.2.1  Independent  Interval  Processing  The  solution  to  the  finite  length  signal  prob¬ 
lem,  extended  to  provide  independent  processing  of  adjacent  intervals,  has  many  advan¬ 
tages  which  render  it  attractive  for  further  study.  Previous  lapped  transform  and  filter 
bank  approaches  require  an  overlap  with  adjacent  data.  One  advantage  of  this  new  ap¬ 
proach  is  that  processing  interval  j  is  not  dependent  on  processing  interval  j  +  1.  If 
real-time  processing  is  required  for  a  given  communication  system,  this  advantage  could 
become  important.  Also,  error  control  is  facilitated  with  this  new  approach,  since  errors  at 
the  boundary  are  not  propogated  into  the  adjacent  interval.  Further,  total  decoupling  of 
adjacent  intervals  allows  for  parallel  processing  of  data.  This  could  also  become  important 
in  real-time  system  applications.  Finally,  processing  M  data  points  with  a  periodically  ex¬ 
tended  signal  requires  only  M  storage  locations  throughout  the  transformation.  Processing 
M  data  points  with  an  overlapped  signal  requires  M  +  /  locations,  where  /  is  the  sum  of 
both  overlaps.  Both  methods  result  in  M  coefficients  in  the  transformation;  however,  the 
previous  overlapping  method  requires  a  larger  buffer.  Further  investigation  of  this  tech¬ 
nique  is  required  to  determine  which  signal  processing  and  communications  applications 
may  benefit  from  its  use. 

5.2.2  Performance  of  the  GLOT  in  the  Presence  of  Noise  The  effects  of  noise  were 
not  considered  in  this  thesis.  In  auy  practical  communication  system,  however,  noise  can 
have  a  significant  impact  on  performance.  Follow-on  efforts  involving  any  application  must 
include  an  analysis  of  the  effects  of  noise.  The  theory  predicts  that  any  noise  could  be 
perfectly  reconstructed,  along  with  the  output  signal,  when  using  all  spectral  coefficients. 
Perfect  reconstruction  of  the  noise  may  not  be  desirable,  however,  and  methods  of  filtering 
and  coding  the  spectral  coefficients  which  eliminate  noise  must  be  investigated.  If  the  noise 
is  somehow  convolved  with  the  signal  of  interest,  it  may  be  liftered  (filtered)  in  the  cepstral 
domain,  since  cepstral  analysis  is  based  on  separation  of  convolved  signals  (6).  In  general, 
the  performance  of  the  generalized  lapped  transform  in  noise  will  depend  on  the  specific 
type  of  noise  present  and  the  specific  application. 
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5.2.3  Cepstral  Processing  The  experiments  discussed  in  Section  4.6.4  indicate  the 
potential  for  a  low  bit  rate  speech  coder  based  on  the  generalized  lapped  transform  and 
cepstral-like  processing.  Advances  in  fast  Fourier  transform  algorithms,  such  as  the  one 
described  in  Chapter  3,  has  made  sophisticated  cepstral  techniques  practical  for  real-time 
applications  such  as  speech  (6:1441).  Additional  research  is  required  to  determine  the  best 
technique  for  compressing  the  ceptral  coefficients,  including  the  number  of  bits  required 
per  coefficient  and  the  amount  of  overhead  resulting  from  position  or  phase  information 
required  for  reconstruction. 

5.2.4  Variable  Orthonormal  Bases  The  theory  provided  by  Suter  and  Oxley  (33) 
allows  for  different  orthonormal  representations  in  adjacent  windows,  even  though  these 
windows  overlap.  This  concept  has  not  been  presented  in  previous  research.  The  pro¬ 
gram  implemented  as  part  of  this  research  did  not  directly  address  multiple  orthonormal 
representations,  as  the  sinusoidal  basis  defined  in  equation  2.3  was  exclusively  used.  Ex¬ 
panding  the  implementation  program  to  encompass  variable  orthonormal  bases,  such  as 
an  orthonormal  polynomial,  must  be  investigated.  A  signal  which  is  not  sinusoidal  in 
some  interval  may  require  many  coefficients  for  representation.  Potentially,  the  number 
of  coefficients  can  be  decreased  if  another  basis  set  more  closely  matched  to  the  signal  is 
chosen  in  that  interval.  For  example,  a  speech  signal  is  generally  sinusoidal  in  nature  dur¬ 
ing  speech  intervals,  but  slowly  varying  in  non-speech  intervals.  Combining  the  features 
of  variable  size  windows  and  variable  orthonormal  bases  could  result  in  an  effective  speech 
compression  algorithm.  A  sinusoidal  basis  would  be  used  to  represent  the  speech  signal  in 
voiced  regions,  while  another  basis  set  could  be  used  to  represent  the  signal  in  unvoiced 
regions.  Although  the  current  implementation  does  not  recognize  these  different  parameter 
intervals,  Switzer  (34)  has  developed  a  program  to  determine  the  voiced/unvoiced  speech 
boundaries.  Combining  the  generalized  lapped  transform  (generalized  to  include  additional 
orthonormal  ba.ses)  and  the  work  of  Switzer  is  an  area  which  merits  further  research. 


Appendix  A.  Ada  Source  Code 


This  appendix  provides  all  source  code  used  during  the  development  and  testing  of 
the  generalized  lapped  transform. 


A.l  Main  Program 


with  Text  Jo; 
with  Complex.Pkg; 
use  Complex.Pkg; 
with  Math-Lib; 
use  Math.Lib; 
with  Vector-Package; 
use  Vector.Package; 
with  Type-Package; 
use  Type-Package; 
with  Print-Package; 
use  Print.Package; 
with  FFT-Pack; 
use  FFT-Pack; 


procedure  GLOT  is 


package  Integer  Jo  is  new  Text  Jo.Integer  Jo(integer); 
package  Float  Jo  is  new  Text  Jo.Float  Jo(float); 

Max-Data-Length  :  constant  :=  100.000; 

Max.Partitions  :  constant  :=  500; 


Big.Daddy. Data.  Vector 
Big-Daddy. Partition.  Vector 

Actual -Data.Length 
Actual -Partitions 

Deriv-Name, 

Alpha-Name, 

Outfile.Name, 

Difference-Name 

Deriv.Name.Length, 
Alpha-Name-Length, 
Outfile-Name.Length. 
Difference.N  ame.Lenglh 

Info-File 


:  Real-Vector  (l..Max.Data.Length); 

:  Partition-Vector  (0. .Max.Partitions); 

:  integer  :=  0; 

:  integer  ;=  0; 


:  string  (I. .30); 


;  integer  range  1..30; 
:  Text  Jo.  File-Type; 
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—  Get.File JData  prompts  the  user  for  the  prepared  info  file  name 

—  and  then  reads  in  the  data  in  the  following  order 

—  (1)  Input  Data  Filename 

—  (2)  Derivative  Output  Filename 

—  (3)  Alpha  Output  Filename 

—  (4)  Output  Data  Filename 

—  (5)  Input  Output  Difference  Filename 

—  (6)  General  Info  Filename 

—  (7)  Number  of  Points  in  the  Data  Vector 

—  (8)  Number  of  required  partitions 

—  (9)  -  (?)  loop  for  input  given  in  (8)  and  get 

—  Boundary  point  then  Overlap  for  each  partition 


procedure  Get.FileJ)ata 


(  Derivative.Filename 
Derivative-Length 
Alpha-Filename 
Alpha-Length 
Output-Filename 
Outname-Length 
Diff.Filename 
Diff-Length 
Infofile 
Partitions 

Number.OfJPartitions 

Data 

Points 


:  in  out  string; 

:  in  out  integer; 

:  in  out  string; 

:  in  out  integer; 

:  in  out  string; 

:  in  out  integer; 

;  in  out  string; 

:  in  out  integer; 

:  in  out  Text  Jo.  File-Type 
:  in  out  Partition.Vector; 

:  in  out  integer; 

:  in  out  Real.Vector; 

:  in  out  integer  )  is 


—  Filenames,  lengths  must  be  in  out  so  they  can 

—  be  written  to  info  file. 

—  Number. Of-Partitions  and  Points  must  be  in  out 

—  Parameters  because  they  are  used  in  a  loop  to 

—  Read  in  the  appropriate  amount  of  data 

Infile, 

Datafile 
TempJnput 
Data-Filename 
Data-Filename-Length 
In-Filename 
In_Filename-Length 
Info.Filename 
Info-Length 

Counter 

begin 

Text  Jo.New-Line; 

TextJo.put(“What  is  the  name  of  the  file  containing  ”); 

TextJo.put(“the  prepared  information  ?  =  ”); 

Text  Jo. Get-Line  (  In.Filename,  In-Filename-Length); 

Text  Jo. Open(  Infile,  TextJo.In.File,  In-Filename(l..In-Filename.Length)); 


;  TextJo.File.Type; 
:  integer; 

:  string(1..30); 

;  integer; 

:  string(1..30); 

:  integer; 

;  string(1..30); 

:  integer; 

;  integer  :=  0; 
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—  Now  I  call  in  the  required  user  data 

Text  Jo. Get-Line  (Infile,  Data.Filename,  Data.Filename-Length); 
Text  Jo. Get-Line  (Infile,  Derivative-Filename,  Derivative-Length); 
Text  Jo. Get-Line  (Infile,  Alpha-Filename,  Alpha-Length); 

Text  Jo.GetXine  (Infile,  Output-Filename,  Outname.Length); 
Text  Jo. Get-Line  (Infile,  Diff-Filename,  Diff-Length); 

Text  Jo. Get-Line  (Infile,  Info-Filename,  Info-Length); 
Integer-Io.get  (Infile, Points); 

Integer  Jo. get  (Infile, Number-Of_Partitions); 


Text  Jo.Greate(Infofile,  Text  Jo. Out-File,  Info_Filename(l..Info-Length)); 

Text  Jo. put  (Infofiie,“This  information  describes  the  outputs”); 

Text  Jo. put  (Infofile,"  produced  from  prepared  data  file  =  ”); 

Text  Jo.  Put-Line  (Infofile,In-Filename(l.. In-Filename-Length)); 

Text  Jo.NewXine  (Infofile); 

TextJo.put  (Infofile, “Input  Data  Filename  >  ”); 

Text  Jo.  Put-Line  (Infofile,  Data-Filename(  1 . .  Data_Filename-Length) ) ; 
Text Jo.NewXine  (Infofile); 

Text  Jo. put(Infofile, “Derivative  Filename  =  ”); 

Text  Jo.  Put-Line(Infofile,  Derivative-Filename  ( 1 . .  Deri  vati  veXength) ) ; 
Text  Jo.NewXine(Infofile); 

Text  Jo. put(Infofile, “Alphas  Filename  =  ”); 

Text  Jo.Put-Line(Infofile,Alpha_Filename(l..AlphaXength)); 

Text  Jo.  NewXine(  Infofile); 

Text  Jo. put(Infofile, “Output  Data  Filename  =  ”); 

Text  Jo.  Put-Line(Infofile,Output-Filename(  1 .  .Outname.Length)) ; 
TextJo.NewXine(  Infofile); 

Text Jo.put(Infofile,“In  Out  Difference  Filename  =  ”); 

Text  Jo.  PutXine(lnfofile,Diff-Filename(l..  Diff-Length)); 

Text  Jo.NewXine(Infofile); 

Text  Jo. put(Infofile,  “Number  of  Data  Points  =  ”); 

Integer-Io.put(lnfofile,  Points,!); 

Text  Jo.NewXine(Infofile);  TextJo.NewXine(Infofile); 

Text  Jo. put(Infofile,  “Number  of  Partitions  =  ”); 
Integer.Io.put(lnfofile,(Number-Of-Partitions  +  1),1); 

Text  Jo.  New_Line(  Infofile); 


for  jin  l..Number_Of-Partitions  loop 

Integer  Jo. get(Infile,Partitions(j).  Boundary); 

Integer  Jo. get(  Infile, Partitions(j). Overlap); 
end  loop; 

Text  Jo.Close(Infile); 

—  now  I  call  in  the  numbers  for  the  vector  Data 

Text  Jo.  Open  (Datafile, Text  Jo.  In.File,  Data_Filename(l..Data_Filename-Length)); 
for  j  in  !.. Points  loop 
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Integer  Jo. Get(Datafile, Temp  Jnput); 

Data(j)  ;=  float(TempJnput);  —  for  integer  inputs 
end  loop; 

Text  Jo.Close(Datafile); 

Partitions(O). Boundary  ;=  1; 

Partitions(O). Overlap  :=  0; 

Partitions(Number-OfJ'artitioiis  +  1). Boundary  :=  Points  +  1; 

Partitions(Number.OUartitions  +  2). Boundary  :=  Points  +  Partitions(l). Boundary; 
Partitions(Numijer.Of Partitions  +  1). Overlap  ;=  0; 

Partitions(Nuinber.OfPartitions  +  2)  Overlap  :=  0; 

Number.OfPartitions  ;=  Number.OfPartitions  +  2; 

—  write  the  first  partition  to  the  end  of  the  data  file 

—  for  wrap  around  processing 
Counter  1; 

for  j  in  (Points  +  1).. (Points  +  Partitions!  1). Boundary)  loop 
Data(j)  :=  Data(Counter); 

Counter  ;=  Counter  +  1; 
end  loop; 

Points  ;=  Points  +  Partitions(l). Boundary; 
er>d  Get.File.Data; 


-  Multiply.By .Window 


procedure  MultiplyPy .Window  (  FullJ)ata.Segment  :  in  Real.Vector; 

FoldedJn.Data  :  in  out  Real.Vector; 

Window.)  :  in  out  Real.Vector)  is 


Pi  :  constant 

EPSILON.)  :  integer 

EPSILON.)_plus.l  ;  integer 

A.)  :  integer 

A.j.plus_l  :  integer 

Window.Boundary .Value  ;  constant 


:=  3.141592654; 

Folded JnJData’first  -  FulLDataPegment’first; 
:=  FullJData.Segment’last  -  FoldedJn.Data’last; 
:=  Folded  Jn.Data’first; 

FoldeaJn.Data’last; 
float  :=  1.0/sqrt(2.0); 


begin 


—  This  is  the  basic  window  -  also  need  complex  window 

- CREATE  WINDOW - 

—  Window  defined  in  Section  2. 1.2.1,  equation  2.8. 


—  left  rising  edge 
if  EPSILON.)  =  0  then 

Window.)(A.))  :=  Window.Boundary. Value; 

else 

for  X  in  (A.)  -  EPSILON  j)..(A.)  +  EPSILON.))  loop 

Window_)(x)  ;=  sin((Pi/(4.0»float(EPSILON.))))  ♦  (float(x-(A.)-EPSILON_))))); 
end  loop; 


end  if; 

—  center 

for  X  in  (A_j  +  EPSILON.j  +  l)..(A-j-plus-l  -  EPSILON.j-plus.l  -  1) 
loop 

Window_j(x)  1.0; 
end  loop; 

—  right  falling  edge 

if  EPSILON.j.plus-l  =  0  then 

Window_j(A_j-plus.l)  :=  Window-Boundary  .Value; 

else 

for  X  in  (A.j.plus.l  -  EPSILON.j.plus.l)..  {A.j-plus.l  +  EPSILON.j.plus.l)  loop 
Window _j(xl  :=  cos((Pi/(4.0*float(EPSILON.j.plus_l)))  * 

(float (x-( A.j.plus.l  -EPSILON.j.plus.l)))); 

end  loop; 
end  if; 


- Multiply  By  Window - 

for  X  in  (A.j)..(A.j  -1-  EPSILON.j)  loop 

Folded Jn.Data(x)  :=  FullJData.Seginent(x)  ♦  Window.j(x)  - 

Full.DataJSegment(2tA.j  -  x)  ♦  Window .j(2*A.j  -  x); 

end  loop; 

for  X  in  (A.j  +  EPSILON.j  -1-  1).. (A.j.plus.l  -  EPSILON.j.plus.l  -  1)  loop 
FoldedJnJData(x)  :=  FullJ)ata.Segment(x)  *  Window .j(x); 
end  loop; 

for  X  in  (A.j.plus.l  -  EPSILON.j.plus.l).. (A.j.plus.l)  loop 

FoldedJn.Data(x)  ;=  FullJ)ata.Segment(x)  »  Window.j(x)  -I- 

FulLData.Segment(2*A.j.plus.l  -  x)  ♦  Window.j(2*A-j.plus.l  -  x) 

end  loop; 


end  Multiply .By.Window; 


—  Compute.Coefficients 

—  This  procedure  takes  in  data  from  0  -  N  or  A.j  -  A.j-(-l  (N-t-1  points)  and 

—  returns  coefficients  stored  from  0  -  N-1  or  A.j  -  (A.j-fl  -  1)  (N  points) 

—  Nth  point  is  used  in  calculations  but  is  not  a  valid  out  point 

—  It  will  be  used  to  store  the  zero’th  coefficient  in  the  next  data  segment 


procedure  Compute.Coefficients  (  Data.S<jgment  :  in  out  Real.Vector; 

Deriv.Output  :  in  out  Text Jo.File.Type)  is 

Two.N  ;  integer  :=  2  *  (Data.Segnient’length  -  1); 

End.FFT.Data  :  integer  :=;  Data.Segment’first  -f-  TwoJV  -  1; 
FFT.Data  :  Complex.Vector(DataT5egment’fir5t.. End.FFT.Data); 
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begin 

—  (1)  Even  Extend  Data-Segment 

—  (See  Vector_Package-Complex-Package  and  Section  2. 1.2. 4,  Step  5). 

FFT.Data  :=  Complex-Of(Even_Extend(Data_Segment)); 

—  (2)  Perform  Inverse  FFT  (See  FFT.Pack  and  Section  2. 1.2. 4,  Step  6). 

FFT  (FFT-Data,  True);  —  includes  1/N  scaling 

for  X  in  Data_Segment’range  loop 

Float  Jo.  put(Deriv_Output,FFT-Data(x).  real); 

TextJo.NewXine(Deriv-Output); 
end  loop; 

—  (3. a)  Assign  zero  coefficient  to  first  point  in  Data_Segment. 

Data-Segment(Data_Segment’first)  :=  FFT-Data(FFT_Data’first).Real  ♦  sqrt(float(Two.N)); 

—  (3.b)  Assign  coefficients  1  to  N-1  back  to  Data.Segment  after  integration 

—  (See  Section  2. 1.2. 4,  Step  7). 

for  k  in  (Data_Segment ’first  +  l)..(Data_Segment’last-l)  loop 

Data-Segment(k)  :=  Data_Segment(k-l)  +  (2.0  *  sqrt(float(Two.N))  *  FFT-Data(k).Real); 
end  loop; 

end  Compute.Coefficients; 


—  Reconstruction  FFT 


procedure  Reconstruction. FFT  (  Jnput.Data  :  in  Real. Vector; 

Output.Data  :  out  Real.Vector  )  is 

:  constant  :=  3.141592654: 

:  float  .=  0.0;  —  float  counter 

:  Complex; 

:  integer  :=  2  *  Input. Data’Iength; 

:  float  :=  float(Two-N); 

:  integer  ;=  Input.Data’first  +  TwoJM  -  I; 

:  CompleX-Vector(Input-Dala’first..End.F'FT.Data); 

begin 


Pi 

X 

Scale.Factor, 
Complex. Factor 

Two.N 
Two.N  J^loat 
End. FFT.Data 
FFT.Data 


—  (1)  Odd  Extend  Input. Data 

—  (See  VectorT’ackage,  Complex-Package  and  Section  2.1.3,  Step  1). 
FFT.Data  :=  CompleX-Of(Odd.Extend(Input.Data)); 

—  (2)  Perform  Inverse  FFT  (See  FFT.Pack  and  Section  2.1.3,  Step  2). 
FFT(FFT.Data,  True); 

X  :=  0.0; 

Scale-Factor. real  :=  sqrt(TwoJv. Float); 

Scale-Factor. imag  ~  0.0; 

for  1  in  Output. Data’range  loop 
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Complex-Factor. real  :=  Cos(Pi*X/Two-N_Float); 
Complex-Factor. imag  :=  Sin(Pi*X/Two-N -Float); 

FFT-Data(l)  :=  Scale-Factor  *  Complex-Factor  ♦  FFT.Data(l); 
X  :=  X  +  1.0;  —  float  counter 
end  loop; 

—  cissign  scaled  imaginary  part  to  Output-Data  1..N 
for  1  in  Output-Data’range  loop 

Output-Data(l)  :=  FFT-Data(l).imag; 
end  loop; 

end  Reconstruction-FFT; 


—  Divide-By.Windows  (See  Section  2.1.3,  Step  3). 


procedure  Divide-By.Windows  (  Data-Segment  :  in  out  Real-Vector; 

Window.Right  :  in  Real-Vector; 

Boundary. Point:  in  integer  )  is 

Counter  :  integer  :=  0; 

Window. Left  :  Real-Vector(Window-Right’range); 

Temp-Data.Segment  :  Real-Vector(Data.Segment’range)  :=  Data.Segment; 

begin 

Window. Left  :=  Reverse-Assignment(Window.Right); 

—  Calculate  Data-Segment(A.j  -  e.j  ..  A.j  -  1) 

Counter  :=  Window.Right’last; 

for  i  in  (Window.Right’first)..(Boundary-Point-l)  loop 

Data-Segment{i)  :=  Window.Left(i)  +  Temp-Data-Segment(i)  - 

Window-Left(Counter)  *  Temp.Data-Segment(Counter); 

Counter  Counter  -  1; 
end  loop; 

—  Calculate  Data.Segment(A.j) 

Data-Segment(Boundary.Point)  :=  Temp.Data-Segment(Boundary-Point)  / 

(2.0  ♦  Window.Left(Boundary-Point)); 

—  Calculate  Data-Segment(A.j  +  1  ..  A.j  +  e.j  -  1) 

Counter  :=  Boundary. Point; 

for  i  in  (Boundary-Point  +  1).. (Window .Right’last  -  1)  loop 
Counter  :=  Counter  -  1; 

Data-Segment(i)  :=  WindoW-Right(Counter)  ♦  Temp.Data-Segment(Counter)  + 
Window .Right(i)  *  Temp-Data-Segment(i); 

end  loop; 

—  Calculate  Data.Segment(A.j  +  e.j  ..  A.j-1  -  e_j-l) 

-  for  i  in  Window. Right ’last. .Temp-Data-Segment ’last  loop 
—  Data-Segment(i)  :=  Temp.Data-Segment(i); 

—  Window  values  here  are  unity  so  no  action  required 
—  If  different  scaling  was  required  this  would  have  to  be 
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—  rewritten  -  best  way  would  be  to  pass  in  another  slice 

—  from  window  and  giving  it  a  new  name  as  input  to  Div_By_Win 
—  end  loop; 

end  Divide-By-Windows; 


—  Do-Work 


procedure  Do-Work  (  Data 

Partitions 

Derivfile-String 

Alphafile-String 

Outfile-String 

DifF-File-String 

Infofile 


in  Real-Vector; 
in  Partition-Vector; 
in  String; 
in  String; 
in  String; 
in  String; 

in  out  TextJo. File-Type)  is 


Pi 


constant  ;=  3.141592654; 


Transformed-Data  ,  T-Data 

Output-Data 

Window 

Median-Data,  Median2-Data 


Real-Vector(Data’range); 

Real.Vector(Data’range); 

Real-Vector(Data’range); 

Real.Vector(Data’range); 


Start-Data,  End.Data, 

Start. Window,  End.Window, 

Last.Output. Point, 

Last-Window-Point, 

Counter, 

Alpha.Compression.Counter  :  integer  :=  0; 

Two-N, 

Temp.Nth.Value  :  float  0.0; 

—  Used  because  Nth  value  needed  in  two  intervals 

—  and  otherwise  would  be  overwritten 
DerivfiIe,Alphafile, 


-  A-j  and  A.j+1 

-  -  A-i  -  e.j  and  A.j+1  +  e.j+1 

-  -  A-j+1  -  e.j+1 

-  -  A-j  +  e.j 


Outfile,  DifF.File 
Alpha.Threshold 
Percent-Remaining 


Text  Jo. File-Type; 
float  ~  0.0; 
float  ;=  0.0; 


begin 

Text.Io. Create  (Derivfile,  Text Jo.Out.File,  Derivfile.String); 
TextJo. Create  (Alphafile,TextJo.Out.File,  Alphafile-String); 
Text  Jo. Create  (Outfile,  Text  Jo.Out.File,  Outfile-String  ); 
Text  Jo. Create  (DifT.File,  TextJo.Out.File,  Diff-file-String); 


Text  Jo.New-Line; 

—  This  threshold  is  optional  depending  on  the  application 
Text  Jo. put(  “What  is  the  desired  coefficient  threshold  ?  ”); 
Float  Jo. get(  Alpha.Threshold); 

Text  Jo.New-Line; 

for  j  in  0..(Partitions’last-l)  loop 

begin  —  each  segment  is  transformed  within  this  loop 
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Assign  j  Dependent  Values- 


Start-Data  :=  Partitions(j). Boundary; 

End-Data  Partitions(j+l). Boundary; 

Start-Window  :=  Partitions{j). Boundary  -  Partitions(j). Overlap; 
End.Window  :=  Partitions(j+l). Boundary  +  Partitions(j+l). Overlap; 

Two-N  :=  float(2  ♦  (End-Data-Start-Data)); 

Last-Output-Point  :=  End-Data  -  Partitions(j+l). Overlap  -  1; 
Last-Window-Point  :=  Start-Data  +  Partitions(j). Overlap; 
Temp.Nth.Value  :=  Transforined-Data(End-Data); 


Text Jo.put(“Working....  Processing  Data  Segment”); 

IntegerJo.put(j+l); 

TextJo.New.Line; 

Multiply-By-Window  (  Data(Start-Window.. End-Window), 

Transformed. Data(Start-Data. .  End.Data) , 
Window(Start.Window..  End-Window)); 

—  Steps  2  and  3  of  Coefficient  Evaluation  (Section  2, 1.2. 4). 

—  Window  added  as  passed  parameter  only  so 

—  it  does  not  have  to  be  recalculated 

—  in  Divide.By -Windows 

Counter  :=  0; 

for  1  in  Start-Data.. End.Data  loop 

Transtormed-Data(l)  :=  Transformed-Data(l)  *  sin((Pi  *  float(Counter))/(Two_N)); 
Counter  :=  Counter  +  1;  —  Counter  goes  from  0  to  N 
end  loop; 

—  Step  4  of  Coefficient  Evaluation  (Section  2. 1.2. 4). 

Compute-Coefficients  (  Transformed.Data(Start-Data.. End-Data), 

Derivfile); 

—  This  routine  performs 

—  (1)  an  even  extension  of  the  data, 

—  (2)  an  inverse  EFT, 

—  (3)  assigns  coefficients  0  -  N-1  to  Translormed.Data 

—  Steps  5,  6,  and  7  of  Coefficient  Evaluation  (Section  2. 1.2.4). 

- Print  Alphas  for  each  segment  to  file - 

for  X  in  Start-Data.. (End.Data-l)  loop 

float  Jo.  put(Alphafile,Transformed-Data(x)); 
text  Jo.newJine(Alphafile); 
end  loop; 


Transformed.Data(End.Data)  :=  Temp-Nth-Value; 
—  Transformed-Data  will  not  be  changed  again  in 
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—  this  loop.  This  assignment  replaces  the  Nth  value 

—  back  to  Transformed-Data  so  it  can  be  used 

—  as  the  Start-Data  point  in  the  following  segment 

Reconstruction-FFT  (  Transformed.Data 

(Start.Data..(  End-Data- 1)), 
Output-Data((Start_Data-t- 1 ) . .  End-Data) ) ; 

—  This  routine  performs 

—  (1)  an  odd  extension  of  the  data, 

—  (2)  an  inverse  FFT, 

—  (3)  assigns  the  scaled  imaginary  result  to 

—  the  Output-Data  values  1  -  N  or  A-j  -I-  1  to  A-j-t-1 

—  This  indicates  that  data  values  0  -  N 

—  create  alphas  from  0  -  N-1  which  in  turn  are  used 

—  to  reconstruct  data  values  from  1  -  N 

—  -  WOW  !!! 

—  Steps  1  and  2  of  Reconstruction  (Section  2.1.3). 

Divide-By. Windows  (  Output-Data(Start-Window..Last-Output-Point), 
Window(Start.Window..  Last-Window -Point), 
Partitions(j) .  Boundary) ; 

if  j=  (Partitions’last  -  1)  then 
begin 

Counter  :=  1; 

for  k  in  Start-Data.. Last.Window.Point  loop 
Output.Data(Counter)  ;=  Output.Data(k); 

Counter  :=  Counter  -I-  1; 
end  loop; 
end;  —  block 
end  if; 

end; - block 

end  loop;  —  j  loop 

TextJo.CIose(Alphafile); 

Text  Jo.  Put.Line(  “Printing  output  to  file...”); 


for  X  in  Partitions(O). Boundary.. (Partitions(Partitions’last  -  1). Boundary  -  1)  loop 
Integer  Jo. put(Outfile,(Integer(Output-Data(x)))); 

Text  Jo.NewXine(Outfile); 

Integer  Jo.put(Difr-File,(Integer(Output.Data(x)  -  Data(x)))); 
TextJo.New_Line(Diff-File); 
end  loop; 

TextJo.New-Line(Infofile); 

TextJo.put(Infofile,“The  threshold  was  ”); 

Float  Jo. put(Infofile,  Alpha-Threshold, 2, 2,0); 

Text  Jo.New-Line(Infofile);TextJo.New_Line(Infofile); 

TextJo.put(Infofile,“No.  of  zeroed  coefficients  was  ”); 

Integer  Jo. put(Infofile,  Alpha-Compression-Counter, 2); 
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Text  jo.NewJLine(Infofile);TextJo.New_Line(Infofile); 

Text  Jo.put(In{ofile,“The  Percent-Remaining  :=  100.0  ♦  (1.0  - 
(float( Alpha-Compression-Counter)  / 
float((Partitions(Partitions’last-l). Boundary  -  1)))); 

Float  Jo. put(Infofile,  Percent-Remaining, 2, 5,0); 

Text  Jo.  N  ew_Line{Infofile); 

Text  Jo.Close(Outfile); 

Text  Jo .  Close(  Diff-F  ile ) ; 

Text  Jo.NewXine; 

Float  Jo. put(  Alpha-Threshold, 3,5,0); 

Text  Jo.NewXine; 

Integer-Io.put( Alpha-Compression-Counter, 3);  Text  Jo.New-Line; 
Float  Jo. put(Percent-Remaining, 3,5,0); 

Text  Jo.  N  ew-Line; 

end  Do-Work; 


begin  —  main 

Get-File.Data  (  Deriv-Name,  Deriv-Name-Length, 

Alpha-Name,  Alpha-Name.Length, 

Outfile.Name,  Outfile-Name-Length, 
DifFerence.Name,  Difference.Name.Length, 

Info.File, 

Big-Daddy-Partition.Vector,  Actual-Partitions, 
Big-Daddy-Data. Vector,  Actual-Data-Length); 

—  By  sending  the  correct  size  arrays  into  Do.Work,  the  array 

—  attributes  can  be  used  to  determine  the  size  rather  than 

—  passing  another  parameter 

Do.Work  (  Big-Daddy. Data. Vector(l.. Actual-Data-Length), 
Big-Daddy-Partition_Vector  (0.. Actual-Partitions), 

Deri  v.N  ame(  1 . .  Deri  v-Name-Length) , 

Alpha.Name(l..  Alpha.N  ame-Length ) , 

Outfile-N  ame(  1 .  .Outfile.Name-Length) , 

Difference.N  ame{  1 . .  Difference.N  ame-Length ) , 

Info-file); 

end;  —  main 
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A .  2  Package  Type-Package 


package  Type.Package  is 

type  Complex  is 
record 

Real  :  float; 

Imag  ;  float; 
end  record; 

type  Partition-Record  is 
record 

Boundary  ;  integer; 

Overlap  :  integer; 
end  record; 

type  Partition-Vector  is  array  (integer  range  <>)  of  Partition-Record; 

—  Partition-Vector  is  an  array  of  Records 
—  Each  record  specifies  the  boundaries 

—  and  the  number  of  points  of  overlap 

type  Direction  is  (Increasing,  Decreasing,  Undefined); 
type  Phase.Type  is  (Pos,  Neg,  Zero); 

type  Phase.Vector  is  array  (integer  range  <>)  of  Phase-Type; 

—  used  only  for  polynomial  fit  compression  technique  outlined  in  Section  4.6.3 

type  Real.Vector  is  array  (integer  range  <>)  of  float; 

type  Complex.Vector  is  array  (integer  range  <>)  of  Complex; 

—  By  using  an  unconstrained  arrays  <>  I  will  be  able  to 

—  pass  array  slices  in  and  out  of  procedures  which  are  sub- 

—  procedures  to  Do-Work.  The  only  arrays  global  to  the 

—  entire  program  are  Big-Daddy-Data- Vector  and 
—  Big-Deiddy-Partition-Vector 

—  All  other  arrays  are  slices  tailored  to  user 

—  input  in  Get-Parameters. 

end  Type-Package; 


A. 3  Package  FFT-Pack 

The  code  for  the  FFT  Package  specification  and  body  are  contained  in  Appendix  A. 


A .  4  Package  Vector-Package 

A.4-1  Vector-Package  Specification 


with  Text  Jo; 
with  Type-package; 
use  Type-Package; 
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with  Complex_Pkg; 
use  Complex.Pkg; 

package  vector-package  is 


function  Odd.Extend(In_Vector  ;  in  Real-Vector)  return  Real-Vector; 
function  Even-Extend(In-Vector  :  in  Real-Vector)  return  Real-Vector; 
function  Reverse-Assignment  (  In-Vector  :  in  Real-Vector)  return  Real-Vector; 
end  vector-package; 


A. 4-2  Vector J^ackage  Body 


package  body  vector.package  is 

package  integerJo  is  new  textJo.integer-io(integer); 

function  Even-Extend(In.Vector  :  in  Real-Vector)  return  Real.Vector  is 
—  See  Section  2. 1.2. 4,  Step  5. 

Two.N  :  integer  :=  2  ♦  (In.Vector’length  -1); 

End.Extended.Data  :  integer  :=  In.Vector’first  +  Two.N  -  1; 
Extended.Data  :  Real-Vector(In.Vector’first.. End.Extended.Data); 

Counter  :  integer  :=  In.Vector’last  -  1; 

begin 

Extended-Data(Extended-Data’first.. In.Vector’last)  :=  In.Vector; 
for  i  in  (In.Vector’last  +  I)..Extended.Data’last  loop 
Extended-Data(i)  :=  In-Vector(Counter); 

Counter  :=  Counter  -  1; 

end  loop; 

return  Extended.Data; 
end  Even.Extend; 


function  Odd.Extend(In-Vector  :  in  Real-Vector)  return  Real-Vector  is 
—  See  Section  2.1.3,  Step  1. 

Two.N  :  integer  :=  2  ♦  In.Vector’length; 

End-Extended-Data  :  integer  ;=  In.Vector’first  +  Two-N  -  1; 
Extended-Data  :  Real-Vector(In-Vector’first.. End-Extended. Data); 

Counter  :  integer  :=  In.Vector’last; 

begin 

Extended-Data(Extended.Data’first..In-Vector’last)  :=  In.Vector; 
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for  i  in  (In.Vector’last  +  l)..Extended_Data’last  loop 
Extended_Data(i)  :=  -  In_Vector(Counter); 
Counter  :=  Counter  -  1; 
end  loop; 

return  Extended-Data; 
end  Odd-Extend; 


function  Reverse_Assignment(In-Vector  :  in  Real.Vector)  return  Real-Vector  is 

Counter  :  integer  :=  In-Vector’first; 

Reversed-Vector  :  Real-Vector(In-Vector’range); 

begin 

for  y  in  reverse  In-Vector’range  loop 

Reversed-Vector(Countet)  In-Vector(y); 

Counter  ;=  Counter  +  1; 
end  loop; 

return  Reversed-Vector; 
end  Reverse-Assignment; 


end  vector-package; 


A .  5  Package  Complex^Pkg 

A. 5.1  Complex.Pkg  Specification 


with  Type.Package; 
use  Type-Package; 

package  Complex-Pkg  is 

—  arithmetic  operations 


function 

function 

function 

function 

function 

function 

function 


(  A,B  ;  Complex)  return  Complex; 

(  A.B  :  Complex)  return  Complex; 

(  A,B  :  Complex)  return  Complex; 

(  R  :  float; 

C  :  Complex)  return  Complex; 

“/"  (  A  :  Complex; 

B  ;  float)  return  Complex; 

“/”  (  A.B  ;  Complex)  return  Complex; 

Negative  (A  ;  Complex)  return  Complex: 


—  conversion  operations 


function  Complex-Of  (R.I  :  float)  return  Complex; 
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function  Complex-Of  (R.Vector  :  Real.Vector)  return  Complex-Vector; 
function  Complex-Of  (R  ;  float)  return  Complex; 
function  Conjugate  (A  :  Complex)  return  Complex; 
procedure  Set.Equal  (  A  :  in  Complex;  B  :  out  Complex); 

—  I/O  operations 

procedure  Get  (A  :  out  Complex); 
procedure  Put  (A  :  in  Complex); 

end  Complex-Pkg; 


A. 5. 2  Complex-Pkg  Body 

with  Math-Lib; 
use  Math-Lib; 
with  Text  Jo; 
use  Text  Jo; 

package  body  Complex-Pkg  is 

package  FltJo  is  new  Float  Jo(float); 
use  FltJo; 

function  “+”  (A,B  :  Complex)  return  Complex  is 
Result  ;  Complex; 
begin 

Result. Real  :=  A. Real  +  B.Real; 

Result. Imag  :=  A.Imag  +  B.Imag; 

Return  Result; 
end  “+”; 

function  (A,B  :  Complex)  return  Complex  is 
Result  :  Complex; 
begin 

Result. Real  :=  A. Real  -  B.Retil; 

Result. Imag  ;=  A.Imag  -  B.Imag; 

Return  Result; 
end 

function  (A,B  :  Complex)  return  Complex  is 
Result  :  Complex; 
begin 

Result. Real  :=  (A.Real*B.Real)  -  (A.Imag*B.Imag); 
Result. Imag  :=  ( A.Real*B.Imag)  -p  (A.Imag+B.Real); 
Return  Result; 
end 

function  (R  :  float;  C  :  Complex)  return  Complex  is 
Result  :  Complex; 
begin 

Result. Real  :=  R  ♦  C.Real; 

Result. Imag  :=  R  *  C.Imag; 

Return  Result; 
end 
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function  “/”  (A  ;  Complex;  B  :  float)  return  Complex  is 
Result  :  Complex; 
begin 

Result. Real  :=  A. Real  /  B; 

Result. Imag  :=  A.Imag  /  B; 

Return  Result; 
end  7”; 

function  “/”  (A,B  :  Complex)  return  Complex  is 
Result  :  Complex; 

Squares  ;  float; 
begin 

Squares  :=  A. Real  **  2  +  B.Imag  ♦♦  2; 

Result. Real  (A. Real  ♦  B.Real  +  A.Imag  ♦  B.Imag)  /  Squares; 
Result. Imag  ;=  (A.Imag  ♦  B.Real  -  A. Real  ♦  B.Imag)  /  Squares; 
Return  Result; 
end  7”; 

function  Negative  (A  :  Complex)  return  Complex  is 
Result  :  Complex; 
begin 

Result. Real  :=  -A. Real; 

Result. Imag  :=  -A.Imag; 

Return  Result; 
end  Negative; 

function  Complex.Of  (R,I  :  float)  return  Complex  is 
Result  ;  Complex; 
begin 

Result. Real  :=  R; 

Result. Imag  :=  I; 

Return  Result; 
end  Complex-Of; 

function  Complex.Of  (R.Vector  ;  Real.Vector)  return  Complex-Vector  is 
Result  :  Complex_Vector(R.Vector’range); 
begin 

for  i  in  R.Vector’range  loop 

Result(i).Real  :=  R.Vector(i); 

Result(i).Imag  :=  0.0; 
end  loop; 

Return  Result; 
end  Complex.Of; 

function  Complex.Of  (R  :  float)  return  Complex  is 
Result  ;  Complex; 
begin 

Result. Real  :=  R; 

Result. Imag  :=  0.0; 

Return  Result; 
end  Complex-Of; 

function  Conjugate  (A  :  Complex)  return  Complex  is 
Result  :  Complex; 
begin 
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Result. Real  :=  A. Real; 

Result. Imag  :=  -  A.Imag; 
return  Result; 
end  Conjugate; 

procedure  Set-Equal  (  A  :  in  Complex;  B  :  out  Complex)  is 
begin 

B.Imag  ;=  A.Imag; 

B.Real  :=  A.Imag; 
end  Set-Equal; 

procedure  Get  (A  :  out  Complex)  is 
begin 

Get(A.Real); 

Get(A.Imag); 
end  Get; 

procedure  Put  (A  :  in  Complex)  is 
begin 

Put  (“(”); 

Put  (A.Real,  1,2,0); 

Put  (V); 

Put  (A. Imag, 5, 2,0); 

Put(T); 

end  Put; 
end  Complex-Pkg; 
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Appendix  B.  Polynomial  Fit  Compression  Technique  Code 


The  following  procedures  were  added  to  support  the  polynomial  fitting  compression 
technique  described  in  Section  4.6.3. 


B.l  Procedure  Maxima-Minima 


procedure  Majdma.Mimma  (  Data 

Threshold 

First_Pos_Maxima, 

First_Neg-Maxima, 

Last_Pos-Maxima, 

Last-Neg.Maxima 

Phase 


;  in  out  Real-Vector; 
:  in  float; 


:  in  out  integer; 

:  out  Phase-Vector)  is 


boolean  :=  true; 
Direction  :=  Undefined; 


integer  :=  0; 


boolean  :=  true; 


First-Pos-Loop,  First-Neg-Loop 
Pos-Slope,  NogJSlope 
Last-Pos-Point.5aved, 

Last-Neg-Point-Saved, 

Last-Pos.Point. Considered, 

Last.N  eg.Point. Considered 
First-Trip.Through-Pos.Decreasing-Slope, 

First. Trip-Through-NegJncreasing-Slope 
begin 

First-Pos. Maxima  :=  0; 

First-Neg-Maxima  :=  0; 

—  These  values  are  set  to  zero  to  be  used  in  this  procedure 

—  to  keep  track  of  the  status  of  First*  and  Last* 

Last-Pos-Maxima  :=  0; 

Last-Neg-Maxima  :=  0; 

—  These  values  are  set  to  zero  for  each  segment  so  they  can  be  checked 

—  upon  exit  from  Maxima-Minima  and  used  as  a  decision  basis  for 

—  reconstructing  using  polynomials 


for  j  in  Data’range  loop 

begin  —  loop 

if  Data(j)  >  Threshold  then 

begin  —  positive  values 

Phase(j)  :=  Pos; 

if  Pos-Slope  =  Increasing  then 

if  Data(j)  >  Data(Last-Pos-Point-Considere  1)  then 
Data(Last.Pos-Point-Considered)  :=  0.0; 
Last-Pos.P'-l.it-Considered  :=  j; 

else 

Last-Pos-Point-Saved  Last-Pos-Point-Considered; 
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Last-Pos-Point-Considered  ;=  j; 

Pos_Slope  :=  Decreasing; 
if  First_Pos_Maxima  /=  0  then 

Last-P'os.Maxima  :=  Last-Pos.Point_Saved; 

—  if  this  point  in  the  algorithm  is  entered  at  least 

—  once  Last.Pos.Maxima  will  change  from  zero  lO  the 
- appropriate  number  and  polynomial  reconstruction 

—  will  take  place  on  the  positive  coefficients. 

—  This  ocnrrs  when  the  slope  changes  from  increasing 

—  to  decreasing, 
end  if; 

end  if; 

elsif  Pos-Slope  =  Decrca'-ing  then 
if  First-Pos-Maxima  =  0  then 

First.Pos.Maxima  ;=  Last_Pos.Point_Saved; 
end  if; 

if  Data(j)  <  Data(Last.Pos.Point-Considered)  then 
Data(Last-Pos-Point.Considered)  :=:  0.0; 
Last.Pos-Point-Considered  :=  j; 
else 

Last-Pos.Point-Saved  :=  Last.Pos-Point.Considered; 
Last.Pos.Point.Considered  :=  j; 

Pos.Slope  :=  Increasing; 
end  if; 

elsif  PosJSlope  =  Undefined  then  —  1st  or  2cd  pos  loop 
if  First.Pos.Loop  then  —  1st  pos  loop 
Last-Pos.Point-Saved  ;=  j; 

First.Pos.Loop  ;=  false; 

^Ise  —  2cd  pos  loop 

if  Data(j)  >  Data(Last.Pos.PointJjaved)  then 
Pos.Slope  :=  Increasing; 
else 

Pos.Slope  :=  Decreasing; 
end  if; 

Last.Pos.Point.Considered  :=  j; 
end  if; 

end  if;  -  Pos  Slope 
end;  —  block  for  positive  values 

elsif  DataC'  <  -Threshold  then 

begin  -  -  negative  values 

Phase(j)  :=  Neg; 

if  Ntg.Slope  =  Increasing  then 
if  First.Neg.Maxima  =  0  then 

First-Neg-Maxima  :=  Last.Ncg.Point.Savcd; 
end  if; 

if  Data{j)  >  Data(Last.Neg_Point.Considered)  then 
Data(Last.Neg.Point.Considerpd)  :=  0.0; 
Last.Neg.Poinl.Considered  :=  j; 
else 
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Last-Neg-PointJSaved  :=  Last-Neg_Point.Considered; 
Last-Neg-Point. Considered  :=  j; 

Neg^lope  ;=  Decreasing; 
end  if; 

elsif  Neg^lope  =  Decreasing  then 

if  Data(j)  <  Data(Last.Neg-Point_Considered)  then 
Data(Last.Neg-Point_Considered)  ;=  0.0; 
Last-Neg-Point_Considered  :=  j; 
else 

Last-Neg_PointJSaved  :=  Last_Neg-Point_Considered; 
Last-Neg-Point_Considered  :=  j; 

Neg-Slope  ;=  Increasing; 
if  First-Neg-Maxima  /=  0  then 

Last-Neg_Maxima  :=  Last.Neg.Point-Saved; 

—  See  comment  above  in  positive  values  algorithm, 
end  if; 
end  if; 

elsif  Neg_Slope  =  Undefined  then  —  1st  or  2cd  neg  loop 
if  First.Neg-Loop  then  —  1st  neg  loop 
Last_Neg-Point-Saved  :=  j; 

First.Neg-Loop  :=  false; 
else  —  2cd  neg  loop 

if  Data(j)  >  Data(Last.Neg_Point_Saved)  then 
Neg..Slope  :=  Increasing; 
else 

Neg-Slope  :=  Decreasing; 
end  if; 

Last-Neg.Point.Considered  :=  j; 
end  if; 

end  if;  —  Neg  Slope 
end;  —  block  for  negative  values 

else 

Data(j)  :=  0.0; 

Phase(j)  :=  Zero; 
end  if; 

end;  —  block 
end  loop; 

end  Maxima.Minima; 


B.2  Procedure  FitSlice 


—  FitSlice 

procedure  FitSlice  (DataJSegment  :  in  out  Real-Vector; 

Phase-Segment  :  in  Phase- Vector; 
fxl -Position  :  in  integer; 

Correct-Phase  :  in  Phase-Type; 

Thresh  :  in  float  )  is 

fxO  :  float  :=  Data-Segment(Data-Segment’first); 

fxl  ;  float  :=  Data-Segment(fxl_Position); 

fx2  :  float  ;=  Data-Segment(Data-Segment’last); 

AlphaO  ;  float  :=  fxO; 

—  Either  fxO  or  AlphaO  could  be  eliminated 
—  but  both  are  provided  for  clarity 
Alphal, 

Alpha2  :  float;  —  To  be  calculated; 

xl  :  float  :=  float(fxl-position  -  Data-Segment’first);  —  xO  =  0 

x2  :  float  :=  float(Data-Segment’last  -  Data-Segment’first); 

X  :  float  :=  1.0;  —  floating  loop  variable 

Factor  :  float  :=  1.0; 

begin 

Alpha2  :=  (  fx2  -  fxO  -  (x2/xl*(fxl-fx0))  )  /  (x2»(x2-xl))  ; 

Alphal  :=  (  fxl  -  fxO  -  (Alpha2*xl*xl)  )  /  xl  ; 

for  j  in  (Data.Segnient’first+l)..(Data_Segment’last-l)  loop 
if  Phase-Segment(j)  =  Correct-Phase  then 

Data-Segment(j)  :=  AlphaO  +  (alphal*x)  +  (Alpha2*x*x)  ; 

—  Check  to  ensure  polynomial  fit  is  not  out  of  allowable  range 
if  Correct-Phase  =  Pos  then 

if  Data-Segment(j)  <  Thresh  then 
Data-Segment(j)  :=  Thresh; 
end  if; 

else  —  Correct-Phase  =  Neg  then 
if  Data-Segment(j)  >  -Thresh  then 
Data-Segment(j)  ;=  -Thresh; 
end  if; 
end  if; 

end  if; 

X  :=  X  -I-  1.0;  —  X  will  always  =  j  -  Data-Segment’first 
end  loop; 

end  Fit-Slice; 
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B.3  Procedure  PolynomiaLFit 


—  Polynomial.Fit 

procedure  PolynomieJ.Fit  (Data  :  in  out  Real-Vector; 

First.Pos-Maxima, 

First-Neg-Maxima, 

Last-Pos-Maxima, 

Last-Neg.Maxima  ;  in  integer; 

Phase  ;  in  Phase.Vector; 

Threshold  :  in  out  float  )  is 

—  Threshold  is  passed  in  only  to  be  used  in  Fit.Slice 

Start  ;  integer  :=  First-Pos.Maxima; 

Middle,  Endslice  ;  integer  ;=  0; 

The_Last_Point.Saved-Was-A.Min  :  boolean  :=  false; 

begin 

if  Last-Pos.Maxima  /=  0  then 

begin  —  positive  reconstruction  is  valid 

for  X  in  (First.Pos-Maxima+l)..Lasl.Pos-Maxinia  loop 
if  Data(x)  >  0.0  then 

—  then  the  point  is  either  a  min  or  a  max 

if  The.Last.Point-Saved-Was.A-Min  then 
Endslice  ;=  x; 

Fit-Slice(Data(Start..  Endslice), 

Phase(Start..  Endslice), 

Middle,  Pos,  Threshold  ); 
The_Last.Poir.t-Saved.Was_A_Min  :=  false; 

Start  :=  Endslice; 

else  —  else  the  last  point 

Middle  :=  x;  —  saved  was  a  max 
The-Last.Point-Savcd-Was-A-Min  :=  true; 

end  if; 

end  if; 

end  loop;  —  positive  reconstruction  loop 

end;  —  pos  recon  valid  block 
end  if; 

—  Now  handle  negative  polynomial  reconstruction 
The.Last-Point-Saved-Was-A-Min  :=  false;  —  re-initialize 
Start  :=  First-Neg-Maxima; 

if  Last-Neg.Maxima  /=  0  then 

begin  —  negative  reconstruction  is  valid 

for  X  in  (First-Neg.Maxima+l)..Last.Neg-Maxima  loop 
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if  Data(x)  <  0.0  then 

—  then  the  point  is  either  a  min  or  a  max 

if  The.Last-Point_Saved_Was_A.^in  then 
Endslice:=  x; 

Fit.Slice(Data(  Start..  Endslice), 

Phase(Start. .  Endslice) , 

Middle,  Neg  ,  Threshold  ); 
The_Last.Point^aved-Was_A._Min  :=  false; 
Start  ;=r  Endslice; 

else  —  else  the  last  point 

Middle  :=  x;  —  saved  was  a  max 
The_Last.Point_Saved.Was_A_Min  :=  true; 

end  if; 

end  if; 

end  loop;  —  negative  reconstruction  loop 

end;  —  neg  recou  valid  block 
end  if; 

end  Polynomial-Fit; 
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Appendix  C.  Independent  Interval  Processing  Code  Modifications 


The  following  procedures  replace  those  given  in  Appendix  A  for  implementation  of 
the  independent  interval  processing  outlined  in  Sections  2.3  and  4.3.6. 


C.l  Get-Fiie.Data 


-  Get_File_Dala  prompts  the  user  for  the  prepared  info  file  name 

-  -  and  then  reads  in  the  data  in  the  following  order 

(1)  Input  Data  Filename 

-  -  (2)  Derivative  Output  Filename 

(3)  Alpha  Output  Filename 

-  -  (4)  Output  Data  Filename 

-  -  (5)  Input  Output  Difference  Filename 

-  -  (6)  General  Info  Filename 

(7)  Number  of  Points  in  the  Data  Vector 
-•  -  (8)  Number  of  required  partitions 

-  (9)  -  (?)  loop  for  input  given  in  (8)  and  get 

-  -  Boundary  point  then  Overlap  for  each  partition 

-  -  All  code  in  Procedure  Gct.File_Data  remains  the  same  except  the  code 
following  the  Datafile  input  loop  which  assignes  the  data  to  vector  Data 

Partitions(O). Boundary  :=  1; 

Partitions(O). Overlap  .=  0; 

Partitions)  Number.OLPartitions  +  I ). Boundary  :=  Points  +  1; 
Partitions) Number.Of.Partition.s  1). Overlap  :=  0; 

Nurnbi  r.Of_Partitions  :=  Nuinber.Of  J’artitions  +  1; 

Points  Points  Partitions) Nuniber.Of-Partitions). Overlap  1; 

end  Get.Filc-Data; 
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C.  2  Divide -By.  Windows 


—  Divide.By. Windows  (See  Section  2.1.3,  Step  3). 


procedure  Divide_By -Windows  (  Data-Segment  ;  in  out  Real.Vector; 

Window  :  in  Real-Vector; 

AJ  :  in  integer  )  is 

Counter  :  integer  :=;  0; 

Temp-Data-Segment  ;  Real-Vector(Data-Scgment’range)  :=  Data_.Segment; 

EJ  :  integer  :=  AJ  -  Data-Segment’first; 

AJplusl  :  integer  ;=  Data-Segment'last  -  EJ; 

begin 

—  (1)  boundary 

Data-Segment( A J)  :=  Temp-Data.Segment(A J)  /  (2.0  ♦  Window(AJ)); 

—  (2)  left  overlap 
Counter  :=  AJ; 

for  i  in  (AJ  +  l)..(A.l  +  EJ  -  1)  loop 
Counter  :=  Counter  -  1; 

Data-Segment(i)  ;=  Window(Countcr)  ♦  Temp-Data-Segment(Counter) 

+  Window(i)  *  Temp.Data-Segment(i); 

end  loop; 

—  -  (3)  middle 

—  for  i  in  ((AJ  +  EJ).. (AJplusl  -  EJ))  loop 

—  llata-Segment(i)  :=  Teinp-Data-Segment(i); 

—  Window  values  here  are  unity  so  no  action  required 
—  If  different  scaling  was  required  this  would  have  to  be 

—  accounted  for 
'  -  end  loop; 

—  -  (4)  right  overlap 
Counter  ;=  Data-Segment’last; 

for  i  in  ((AJplusl  -  EJ  +  1).. (AJplusl  -  I))  loop 
Counter  :=  Counter  -1; 

Data-Segment(i)  :=  Window(i)  ♦  Temp-Data-Segment(i)  - 

Window(Counter)*Temp-Data-Seginent(Counter); 

end  loop; 

end  Divide-By-Windows; 
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C.3  Do.Work 


—  Do-Work 


procedure  Do-Work  (  Data 

Partitions 

DerivfileJString 

Alphafile-String 

Outfile-String 

Diff-File-String 

Infofile 

Pi 

Transformed-Data 

Window 

aj,  aj.plusl,  ej, 
ajJess.ej,  aj.pluslJess.ej, 
aj.plus.ej,  aj.plusl-plus.ej, 
Counter, 

Alpha.Compression.Counter 

Two.N, 

Temp-N  th.Value 

Derivfile,Alphafile, 

Outfile,  Diff.File 
Two.N, 

Alpha.Threshold, 

Percent. Reiiiainiiig 


:  in  Real  .Vector; 

:  in  Partition-Vector; 

:  in  String; 

:  in  String; 

:  in  String; 

:  in  String; 

:  in  out  Text  Jo.File.Type)  is 

:  constant  :=  3.141592654; 

:  Rea!-Vector(Data’range); 

:  Real_Vector(Data’range); 


:  integer  ;=  0; 
:  float  :=  0.0; 


:  TextJo.File.Type; 


:  float  .=  0.?; 


begin 

Text  Jo. Create  (Derivfile,  TextJo.Out.File,  Derivfile.String); 
Text  Jo. Create  (Alphafile,  TextJo.Out.File,  AlphafileJString); 
Text  Jo. Create  (Outfile,  TextJo.Out.File,  Outfile.String  ); 
TextJo. Create  (Diff.File,  TextJo.Out.File,  DiffJile.String); 

Text  Jo.NewXine; 

Text  Jo. put(  “What  is  the  desired  absolute  threshold  ?  ”); 
Float  Jo. get(  Alpha-Threshold); 

Text  Jo.NewJjine; 


for  j  in  0..(Partitions’last-l)  loop 

begin  —  each  segment  is  transformed  within  this  loop 

- Assign  j  Dependent  Values - 

ej  :=  Partitions(j). Overlap; 
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aj 

aj-plus.ej 
aj  Jess_ej 


=  Partitions(j). Boundary; 
=  aj  +  ej; 

=  aj  -  ej; 


aj-plusl 

aj.pl  us  IJess-ej 

aj-plusl.plus.ej 


;=  Partitions(j4-l). Boundary  ; 
=  aj.plusl  -  ej; 

=  aj.plusl  +  ej; 


TwoJ^ 


:=  float(2  *  (aj.plusl-aj)); 


Text  Jo. put(  “Working....  Processing  Data  Segment”); 

Integer  Jo. put(j+l ); 

Text  Jo.  New. Line; 

Transfotmed.Data(ajJess_ej..aj-l)  ;=  Dala( aj.plusl  Jess-ej..aj.plusl-l); 
Transformed-Data(aj..aj.plusl-l)  ;=  Data{aj..aj_plusl-1); 

Transformed. Data( aj.plusl.. aj.plusl. plus.ej)  :=  Data(aj..aj.plus.ej); 

Multi  ply. By  .Window  (  Transformed.Data(ajJess.ej..  aj.plusl  .plus.ej), 

Transformed  .D  at  a(  aj . .  aj  .pi  us  1 ) , 

Window(aj  Jess  .ej.. aj.plusl. plus.ej)); 

—  Steps  2  and  3  of  Coefficient  Evaluation  (Section  2. 1.2. 4). 

—  Window  added  as  passed  parameter  only  so 

—  it  does  not  have  to  be  recalculated 

—  in  Divide.By .Windows 

Counter  :=  0; 

for  1  in  aj. .aj.plusl  loop 

Transformed. Data(l)  :=  Transformed-Data(l)*sin((Pi*float(Counter))/(Two.N)); 
Counter  :=  Counter  +1;  —  Counter  goes  from  0  to  N 
end  loop; 

—  Step  4  of  Coefficient  Evaluation  (Section  2. 1.2. 4). 

Compute.Coefficients  (  Transformcd.Data(aj.. aj.plusl), 

Derivfile); 

—  This  routine  performs 

—  (1)  an  even  extension  of  the  data, 

—  (2)  an  inverse  FFT, 

—  (3)  assigns  coefficients  0  -  N-1  to  Transformed.Data 

—  Steps  5,  6,  and  7  of  Coefficient  Evaluation  (Section  2. 1.2. 4). 

- Print  Alphas  for  each  segment  to  file - 

for  X  in  aj..(aj.plHsl-l)  loop 

if  ab8(Transformcd.Data(x))  <  Alpha.Threshold  then 
Transformed-nata(x)  ;=  0.0; 

Alpha.Compression.Counter  :=  Alpha.Compression.Counter+1; 
end  if; 

Float  Jo.  put  (Alphafile,Transformed.Data(x)); 

Text  Jo.New  Jine(  Alphafile); 


end  loop; 


Reconstruction-FFT  (  Ttansfotmed_Data(  aj..(aj_plusl-l)), 

Transformed-Data(  ( aj+ 1 ) . . aj -plusl )  ) ; 

Transforined_Data(aj Jess_ej..aj)  :=  Transforined_Data(aj. plusl -less.ej.  aj-plusl); 

Transforined_Data(aj.plusl  +  l..aj-plusl-plus_ej)  :=  Transformed-Data(aj  +  l.,aj_plus_ej); 

Divide.By-Windows  (  Transfornied_Data(aj Jess_ej.,aj_plusl-plus_ej), 

Window  (aj  Jess.ej . .  aj  .pi  us  1  _pl  us_ej ) , 

Parti  tions(j). Boundary); 


for  X  in  (aj..(aj.plusl-l))  loop 

Integer.Io.put(Outfile,(Integer(Transformed_Data(x)))); 

Text_Io.NewJine(Outfile); 

Integer.Io.put  (DifF-File,(Integer(Transformed-Data(x)  -  Data(x)))); 
Text_Io.New-Line(DifF.File); 
end  loop; 
end;  —  block 
end  loop; 

Text  Jo.Close(Alphafile); 

Text  Jo,Close(Outfile); 

Text  Jo,Close(DifF.File); 

TextJo.New_Line(Infofile); 

Textio.put(Infofile,“The  threshold  was  >  ”): 
FloatJo.put(Infofile,Alpha.Threshold,2,2,0); 

Text  Jo.NewJLine(Infofde);Text  jo.New_Line(Infofile); 

TextJo.put(Infofile, “Number  of  zeroed  coefficients  was  >  ”); 

Integer  Jo.  put  (Infofile,Alpha_Compression-Counter,2); 

Text  Jo.New_Line(Infofile);Text  Jo.New_Line(Infofile); 

TextJo.put(lnfofile,“The  %  compression  was  >’^); 

Percent.Remaining  :=  100.0  ♦  (1.0-  (float(Alpha.Compression_Counter)  / 
float((Partitions(Partitions’last-l). Boundary  -  1)))); 

Float  Jo. put(Infofile,  Percent.Remaining,2,5,0); 

TextJo.NewXine(Infofile); 

Text  Jo.NewJjine; 

F'oat  Jo.  put  (Alpha-Threshold, 3, 5,0); 

'I'ext  Jo.NewJ.ine; 

lnteger-Io.put(Alpha.Compression.Coiinter,3); 

Text  Jo.New_Line; 

Float  Jo. put(Percent.Remaining,3, 5,0); 

Text  Jo.  New -Line; 

end  Do-Work; 


C-5 


Vita 


Captain  Brian  Dean  Raduenz  was  born  in  Minneapolis,  Minnesota  on  21  March 
1966.  He  graduated  from  Eveleth  High  School,  Eveleth,  MN  in  1984,  and  entered  the 
United  States  Air  Force  Academy.  He  graduated  from  the  Air  Force  Academy  in  1988 
with  a  Bachelor  of  Science  Degree  in  Electrical  Engineering,  and  a  commission  as  a  Second 
Lieutenant  in  the  U.S.  Air  Force.  Capt  Raduenz  was  then  assigned  to  Hanscom  Air  Force 
Base,  Bedford,  MA.  where  he  served  as  a  project  manager  on  the  Modular  Control  System 
program  (ESD/TC)  until  April,  1991.  He  entered  the  Air  Force  Institute  of  Technology 
in  .May.  1991.  Upon  graduation,  Capt  Raduenz  will  be  assigned  to  Wright- Patterson  Air 
Force  Base.  .\SC/EN. 


Permanent  address:  3306  Cedar  Island  Drive 
Eveleth,  Minnesota  55734 


VITA-1 


Bibliography 


1.  Akansu,  Ali  N.  and  Frank  E.  Wadas.  “On  Lapped  Orthogonal  Transforms,”  IEEE 
Transactions  on  Signal  Processing,  SP-^<?(2):439-443  (February  1992). 

2.  Auscher,  Pascal,  et  al.  “Local  Sine  and  Cosine  Bases  of  Coiffman  and  Meyer.” 
Wavelets  -  A  Tutorial  in  Theory  and  Applications  edited  by  C.  K.  Chui,  237-256, 
Academic  Press,  1992. 

3.  Beylkin,  G.,  et  al.  “Fast  Wavelet  Transforms  and  Numerical  Algorithms  I,”  Commu¬ 
nications  on  Pure  and  Applied  Mathematics,  A'^L/V;141-183  (1991). 

4.  Beylkin,  G.,  et  al.  “Wavelets  in  Numerical  Analysis.”  Wavelets  and  Their  Applications 
edited  by  M.  B.  Ruskai,  et  al.,  181-210,  Boston,  MA:  Jones  and  Bartlett,  1992. 

5.  Bradley,  Jonathan  N.,  et  al.  “Reflected  Boundary  Conditions  for  Multirate  Filter 
Banks.”  Proceedings  of  the  lEEE-SP  International  Symposium  on  Time- Frequency 
and  Time-Scale  Analysis.  307-310.  4-6  October  1992. 

6.  Childers,  Donald  G.,  et  al.  “The  Cepstrum:  A  Guide  to  Processing.”  Proceedings  of 
the  IEEE.  1428-1443.  October  1977. 

7.  Coifman,  Ronald,  “Test  of  the  Bell.”  Private  Communication,  April  1992. 

8.  Coifman,  Ronald  et  Meyer,  Yves.  “Remarques  sur  I’analyse  de  Fourier  a  fenetre,” 
C.R.  Academie  Sci.  Paris,  t.312{SeT\e  I):259-261  (1991). 

9.  Daubechies,  Ingrid.  “Orthonormal  Bases  of  Compactly  Supported  Wavelets,”  Comm. 
Pure  and  Applied  Math,  .^1:909-966  (1988). 

10.  de  Queiroz,  Ricardo  L.  “Subband  Processing  of  Finite  Length  Signals.”  Proceedings  of 
International  Conf.  on  Accoustics,  Speech,  and  Signal  Processing.  lV-613  -  IV-616. 
23-26  March  1992. 

11.  Dick,  David.  “Despite  Resistence,  Ada  Gains  Respectability,”  Signal  Magazine, 
^5(ll):92-94  (1991). 

12.  Fergu.son,  Warren  E.  Jr.  “A  Simple  Derivation  of  Glassman’s  General  N  Fast  Fourier 
Transform,”  Computers  and  Mathematic.s  with  Applications,  8(6):401-411  (1982). 

13.  Glassman,  J.  A.  “A  Generalization  of  the  Fast  Fourier  Transform,”  IEEE  Trans. 
Comput.,  C-/P:105-116  (I970p 

14.  Harris,  Fredric  J.  “On  the  Use  of  Windows  for  Harmonic  Analysis  with  the  Discrete 
Fourier  Transform.”  Proceedings  of  the  IEEE,  66.  51-83.  January  1978. 

15.  Kabrisky,  Matthew.  Thesis  Discussions,  August  1992. 

16.  Karlsson,  Gunnar  and  Martin  Vetterli.  “Extension  of  Finite  Length  Signals  for  Sub- 
Band  Coding,”  Signal  Processing,  7(2);161-168  ( June  1989). 

17.  Mallat,  Stephane  and  Zhifeng  Zhang.  Matching  Pursuits  with  Time-Frequency  Dictio¬ 
naries.  Technical  Report  619,  Courant  Institute  of  Mathematical  Sciences.  October 
1992. 


BIB-1 


18.  Malvar,  Henrique  S.  “Lapped  Transforms  for  Efficient  Transform/Subband  Coding,” 
lEEFy  Transactions  on  Accoustics,  Speech,  and  Signal  Processing,  58(6):969 -978  (June 
1990). 

19.  Malvar,  Henrique  S.  “Extended  Lapped  Transforms:  Properties,  Applications, 
and  Fast  Algorithms,”  IEEE  Transactions  on  Signal  Processing,  .^0(11):270.3-2714 
(November  1992). 

20.  Malvar.  Henrique  S.  Signal  Processing  with  Lapped  Transforms.  Boston  :  Artecli 
House,  1992. 

21.  Malvar,  Henrique  S.  and  David  H.  Staelin.  “The  LOT:  Transform  Coding  Without 
Blocking  Effects,”  IEEE  Transactions  on  Accoustics,  Speech,  and  Signal  Proces.sing. 
.77(4):.5.33-559  (April  1989). 

22.  NIST.  The  D.ARP.\  TIMIT  Acoustic-Phonetic  Continuous  Speech  Corpus  (TIMTT): 
Training  and  Test  Data  and  Speech  Header  .Software,  October  1990. 

23.  Nuri,  Veyis  and  Roberto  H.  Bamberger.  “A  Theoretical  Framework  for  the  Analysis 
and  Design  of  Size- Limited  Multirate  Filter  Banks.”  Proceedings  of  the  lEEE-SP  In- 
ternatiorml  Symposium  on  Time-Frcejuency  and  Time-Scale  Analysis.  311-314.  4-6 
October  1992. 

24.  Oktay,  Alkin.  PC-DSP.  New  Jersey:  Prentice  Hall,  1990.  3|”  IBM  Version. 

2.').  Oppenheim,  Alan  V.  and  Ronald  W.  Schafer.  Discrete  Time  Signal  Processing.  New 
Jersey:  Prentice  Hall,  1989. 

26.  Parsons.  Thomas  W.  Voice  and  Speech  Processing.  McGraw  Hill,  1987. 

27.  Princen,  John  P.  and  Alan  Bernard  Bradley.  “Analysis/Synthesis  Filter  Bank  Design 

Based  on  Time  Domain  Aliasing  Cancellation,”  IEEE  Transactions  on  Accoustics, 
Speech,  and  Signal  Processing,  (5):  11.53-1 161  (October  1986). 

28.  Quackenbush,  Schuyler  R.,  et  al.  Objective  Measures  for  Speech  Processing.  New 
Jersey:  Prentice  Hall,  1988. 

29.  Raduenz,  Brian  D.,  et  al.  “Analysis  of  an  Ada  Ba.sed  Version  of  Glassman's  General 
N  Point  Fast  Fourier  Transform,”  Computers  and  Mathematics  with  Applications 
(Submitted  July  1992,  Accepted  October  1992). 

.30.  Singleton,  R.  C.  “An  Algorithm  for  Computing  the  Mixed  Radix  Fast  Fourier  Trans¬ 
form.”  IEEE  Trans.  Audio  Electroacoust.,  A 7^-/7:93-103  (1969). 

31.  StremI'  ' .  Ferrel  G.  Introduction  to  Communication  Systems  (2nd  Edition).  Addi.son- 
Wesley  Publishing  Co.,  1982. 

32.  Suter,  Bruce  W.  Thesis  Discussions,  January-August  1992. 

33.  Suter.  Bruce  W.  and  Mark  E.  Oxley.  “On  Variable  Overlapped  Windows  and 
Weighted  Orthonormal  Bases,”  IEEE  Tramsactions  on  Signal  Processing  (Subndtted 
April  19.j2). 

34.  Switzer,  Shane.  Frequency  Domain  Speech  Coding.  MS  thesis,  AFIT/GE/ENG/91D, 
School  of  Engineering,  Air  Force  Institute  of  Technology  (AU),  Wright-Patter.son  AFB 
OH,  December  1991. 


BIB  2 


3  j.  Vaidyanathan,  P.  P.  Multirate  Systems  and  Filter  Banks.  New  Jersey:  Prentice  Hal!, 
October  1992. 

36.  Van  Loan,  Charles.  Computational  Frameworks  for  the  Fast  Fourier  Transform. 
Philidelphia:  Society  for  Industrial  and  Applied  Mathematics,  1992. 

37.  Vetterli,  Martin  and  Didier  Le  Gall.  “Perfect  Reconstruction  FIR  Filter  Banks:  Some 
Properties  and  Factorizations,”  IEEE  Transactions  on  Accoustics,  Speech,  and  Signal 
Processing.,  57(7):1057-1071  (July  1989). 

38.  Weinstein,  Clifford  J.  “Opportunities  for  Advanced  Speech  Processing  in  Military 
Computer-Based  Systems,”  Proceedings  of  the  IEEE,  7P(11):1627-1639  (November 
1991). 


BIB-3 


December  1992 


Masters’  Thesis 


DIGITAL  SIGNAL  PROCESSING  USING  LAPPED  TRANSFORMS 
WITH  VARIABLE  PARAMETER  WINDOWS 
AND  ORTHONORMAL  BASES 


Brian  D.  Raduenz 


Air  Force  Institute  of  Technology,  WPAFB  OH  45433-6583 


AFIT/GE/ENG/92D-30 


Jon  Sjogren,  Ph.D. 
AFOSR/NM 

Bolling  AFB,  DC  20332-6448 


Distribution  Unlimited 


This  thesis  develops  and  evaluates  a  number  of  new  concepts  and  tools  for  the  analysis  of  signals  using  variable 
overlapped  windows  and  orthonormal  bases.  Windowing,  often  employed  as  a  spectral  estimation  technique,  can 
result  in  irreparable  distortions  in  the  transformed  signal.  By  placing  conditions  on  the  window  and  incorporating 
it  into  the  orthonormal  representation,  any  signal  distortion  resulting  from  the  transformation  can  be  eliminated 
or  cancelled  in  reconstruction.  This  concept  is  critical  to  the  theory  underpinning  this  thesis.  As  part  of  this 
evaluation,  a  tensor  product  based  general  IV-point  fast  Fourier  transform  algorithm  was  implemented  in  the 
DOD  standard  language,  Ada.  The  most  prevalent  criticism  of  Ada  is  slow  execution  time.  This  code  is  shown  to 
be  comparable  in  execution  time  performance  to  the  corresponding  FORTRAN  code.  Also,  as  part  of  this  thesis, 
a  new  paradigm  is  presented  for  solving  th"  finite  length  data  prob’-'m  Eissociated  with  filter  banks  and  lapped 
transforms.  This  result  could  have  significant  importance  in  manv  Air  Force  applications,  such  as  processing 
images  in  which  the  objects  of  interest  are  near  the  borders.  Additionally,  a  limited  number  of  experiments  were 
performed  with  the  coding  of  speech.  The  results  indicate  the  lapped  transform  evaluated  in  this  thesis  has 
potential  as  a  low  bit  rate  speech  coder. 
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